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1.0  SUMMARY 


The  objective  of  this  effort  is  to  facilitate  reliable  assessment  of  the  durability  and  damage 
tolerance  of  composite  aircraft  in  accordance  with  the  strategic  needs  of  the  Air  Force.  This 
report  outlines  three  critical  tasks  concerned  with  the  development  and  implementation  of  an 
advanced  and  novel  modeling  framework  for  predicting  damage  accumulation  in  composite 
structures  subjected  to  monotonic  and  fatigue  loadings.  The  modeling  framework  accounts  for  the 
physical  mechanisms  of  composite  failure,  fatigue  failure  in  particular,  and  has  been  validated 
based  on  experimental  findings.  The  paragraphs  below  describe  each  task  and  provide  a  guide  to 
the  detailed  information  presented  in  this  report. 

In  Task  1,  a  reduced  order  computational  homogenization  method  to  accurately  capture  the 
physics  of  failure  in  heterogeneous  materials  is  developed  (Section  3.2).  This  multiscale  method 
provides  a  computationally  efficient  modeling  framework  for  predicting  failure  of  composite 
materials  by  directly  accounting  for  damage  accumulation  in  the  constituent  phases  and  in  the 
interface  between  constituents.  An  example  problem  is  included  to  demonstrate  the  capabilities  of 
this  technique  in  the  context  of  particle  reinforced  composites. 

In  Task  2,  a  new  multiple  spatio-temporal  scale  modeling  technique  for  damage 
accumulation  in  composites  subjected  to  cyclic  loading  (Section  3.3)  is  developed.  In  addition  to 
capturing  the  multiple  spatial  scale  effects  addressed  in  Section  3.2,  this  model  addresses  the  scale 
disparity  between  the  characteristic  time  of  a  single  loading  cycle  and  the  total  lifetime  of  a 
composite  structure.  The  presented  model  offers  accurate  life  prediction  while  only  requiring  the 
explicit  simulation  of  a  small  subset  of  loading  cycles  over  a  the  life  of  a  structure. 

Task  3  consisted  of  validation  of  the  model  predictions  against  experimental  data  and 
gaining  an  understanding  of  the  evolution  of  damage  processes  in  carbon  fiber  reinforced 
polymers  (CFRP)  under  monotonic  and  cyclic  loadings  (Section  4).  The  experiments  explored 
failure  mechanisms  using  nondestructive  damage  inspection  techniques  and  provided 
experimental  data  for  both  calibration  of  all  model  parameters  and  validation  of  the  calibrated 
model  response.  The  combined  experimental  and  computational  study  illuminated  the  modes  and 
progression  of  failure  in  monotonically  and  cyclically  loaded  CFRP  composites. 
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2.0  INTRODUCTION 


There  is  little  disagreement  that  modeling  failure  in  brittle  composites  is  of  great 
importance  to  fully  realizing  the  potential  gains  such  materials  offer.  Brittle  composites,  such  as 
carbon  fiber  reinforced  polymers,  often  possess  beneficial  properties  including  a  high  strength  to 
weight  ratio  and  high  fatigue  durability.  These  properties  lend  to  the  usage  of  composites  in 
modem,  high-performance  structures.  However,  without  a  capability  to  model  and  predict  failure, 
large  factors  of  safety  must  be  incorporated  into  designs  preventing  optimum  utilization  of  these 
materials.  This  leads  to  a  high  demand  for  a  failure  modeling  capability  especially  within  the 
aerospace  industry  where  perfonnance  gains  are  paramount. 

The  creation  of  a  predictive  failure  model  for  composite  structures  presents  several 
challenges.  First,  a  wide  variety  of  failure  mechanisms  can  occur  within  the  microstructure  of  a 
composite  material.  These  include  diffuse  microcracking  within  the  matrix,  fiber/matrix 
debonding,  delamination,  fiber  kinking,  fiber  buckling,  and  fiber  fracture.  This  multiplicity  of 
failure  mechanisms  interact  contributing  to  the  ultimate  global  failure  of  a  composite  structure. 
Second,  a  size  disparity  exists  between  the  size  scale  of  constituent  materials  where  failure  initiates 
and  grows  and  the  total  size  of  the  composite  structure.  The  small  scale  geometry  of  the 
intenningled  constituent  materials  cannot  be  resolved  for  an  entire  composite  structure  as  this 
would  require  excessive  amounts  of  computer  memory  and  computational  power.  A  predictive 
model  for  failure  in  composite  structures  must  address  each  of  these  challenging  aspects. 

One  of  the  main  foci  of  this  report  is  predictive  modeling  of  composite  materials  when 
subjected  to  fatigue  loading.  Composites  are  well-known  to  have  high  fatigue  durability,  but  they 
do  fail  in  fatigue.  As  such,  designers  must  account  for  fatigue  failure  which  can  be  both  sudden  and 
catastrophic.  A  predictive  model  for  fatigue  failures  may  aid  designers  to  realize  the  full 
performance  potential  offered  by  composites.  Modeling  fatigue  failure  in  composites  comes  with 
an  additional  challenge.  Much  like  the  previously  discussed  spatial  scale  disparity,  the  duration  of 
a  single  loading  cycle  is  often  significantly  shorter  than  the  total  lifetime  of  a  composite  structure. 
If  millions  of  loading  cycles  are  required  to  induce  global  failure,  explicitly  modeling  each  cycle  of 
loading  would  be  computationally  intractable.  To  model  fatigue  failure  in  composites,  the 
computational  modeling  community  must  treat  the  temporal  scale  disparity  alongside  the  other 
challenges  of  composite  failure  modeling  mentioned  above. 

The  field  of  multiscale  modeling  addresses  each  of  the  stated  challenges  including  both  the 
multiplicity  of  failure  modes  and  the  scale  disparities  in  space  and  time.  Multiscale  modeling  takes 
as  its  aim  the  development  of  methodologies  that  link  phenomena  occurring  at  different  scales. 
The  previously  discussed  scale  disparity  between  a  composite  structure  and  the  small  scale 
interplay  of  the  composite's  constituents  provides  an  example  of  a  problem  with  multiple  spatial 
scales.  An  effective  multiscale  method  for  modeling  failure  in  composites  should  link  the  global 
structural  response  to  the  multiplicity  of  microscale  composite  failure  mechanisms  in  a 
computationally  tractable  manner.  If  fatigue  failure  is  considered,  the  scale  disparity  between  a 
single  loading  cycle  and  the  lifetime  of  a  structure  provides  an  example  of  a  problem  with  multiple 
temporal  scales.  A  tractable  multiscale  model  for  fatigue  failure  in  composites  must  allow 
prediction  of  a  composite  structure's  fatigue  life  without  requiring  the  resolution  of  millions  of 
loading  cycles.  If  such  a  multiscale  approach  can  be  developed,  multiscale  modeling  provides  a 
promising  research  path  directly  attacking  the  most  difficult  aspects  of  modeling  failure  in 
composites. 
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This  report  presents  multiple  spatial  and  temporal  scale  methodologies  that  address  the 
problem  of  modeling  failure  in  composites  subjected  monotonic  and  fatigue  loadings.  The  spatial 
and  temporal  multiscale  methodologies  are  designed  to  operate  simultaneously.  Within  the 
multiple  spatial  scale  methodology,  various  failure  modes  can  be  naturally  incorporated  including 
fiber  failure,  fiber/matrix  debonding,  matrix  cracking,  and  delamination.  To  obtain  computational 
tractability,  a  spatial  order  reduction  and  an  adaptive  time  stepping  technique  are  devised.  These 
reductions  introduce  some  error,  but  high  computational  efficiency  is  gained  while  still 
maintaining  the  important  features  of  the  problem.  These  new  multiscale  methodologies  attempt  to 
directly  address  the  challenges  of  modeling  failure  in  brittle  composites. 

The  remainder  of  this  report  is  organized  as  follows:  Section  3  presents  the  methods 
developed  in  pursuit  of  our  goals  along  with  assumptions  and  experimental  procedures.  Section  4 
presents  the  results  and  discussion  of  an  intensive  computational  and  experimental  study 
conducted  on  the  CFRP  composite,  IM7/977-3.  The  study  was  conducted  to  better  understand 
damage  growth  in  composites  subjected  to  mono  tonic  and  fatigue  loadings  and  to  validate  the 
proposed  modeling  techniques.  Section  5  presents  the  conclusions  of  this  work. 
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3.0  METHODS,  ASSUMPTIONS,  PROCEDURES 

This  section  presents  all  of  the  experimental  procedures  and  numerical  methods  used  to 
produce  the  results  and  discussion  of  Section  4.  We  state  all  assumptions  alongside  each  of  these 
procedures  and  methods.  In  this  work,  our  aim  is  to  produce  accurate  and  computationally  efficient 
multiscale  models  for  damage  accumulation  in  CFRP  composites  undergoing  both  monotonic  and 
fatigue  loadings.  A  multiple  spatio-temporal  modeling  framework  has  been  developed  in 
accordance  with  this  aim  and  will  be  presented  in  this  section.  In  addition,  an  extensive 
experimental  program  was  undertaken  to  both  validate  the  modeling  framework  and  to 
experimentally  ascertain  the  nature  of  damage  growth  in  composites  using  non-destructive  testing. 
Section  3.1  states  the  experimental  procedures  employed  in  the  experimental  program.  Section  3.2 
presents  a  multiple  spatial  scale  technique  for  modeling  damage  growth  in  monotonically  loaded 
composites.  In  Section  3.3,  we  extend  the  method  of  Section  3.2  by  incorporating  multiple 
temporal  scales  allowing  efficient  modeling  of  damage  growth  in  composites  subjected  to  fatigue 
loadings.  Both  Sections  3.2  and  3.3  provide  verification  studies  for  the  developed  techniques. 

3.1  Experimental  Procedures 

A  series  of  monotonic  and  fatigue  tension  tests  were  conducted  on  the  graphite  fiber 
reinforced  epoxy,  IM7/977-3.  In-situ  acoustic  emission  monitoring  was  conducted  in  order  to 
characterize  damage  propagation  with  increasing  load.  X-ray  radiography  and  X-ray  computed 
tomography  were  used  periodically  to  visually  inspect  the  type,  location,  and  extent  of  internal 
damage. 

3.1.1  Material  Fabrication 

Quasi-isotropic  panels  were  hand  laid  from  unidirectional  preimpregnated  IM7/977-3 
graphite  epoxy.  They  were  cured  in  an  autoclave  at  a  temperature  of  1 77°  C  and  a  pressure  of  689 
kPa.  After  cure,  the  panels  were  cut  into  multiple  test  specimens  with  nominal  dimensions  of  25.4 
mm  x  2  mm  x  254  mm.  The  specimens  consisted  of  16  plies  with  [+45,0,-45,90]2s  layup.  The 

mean  and  standard  deviation  of  the  fiber  volume  fraction  were  detennined  to  be  66.6%  and  2.5%, 
respectively,  by  acid  digestion  testing. 

3.1.2  Testing 

Two  sets  of  mono  tonic  tension  tests  were  conducted  on  an  MTS  universal  testing  machine 
according  to  ASTM  D3039  [3].  The  first  set  was  conducted  at  a  constant  displacement  rate  of  1.27 
mm/min  to  obtain  the  average  mechanical  properties  of  0“  and  90°  unidirectional  composite 
specimens.  The  second  set  of  tests  was  conducted  on  a  quasi-isotropic  specimen  in  order  to 
thoroughly  characterize  the  quantity  and  location  of  damage  progression  as  a  function  of  load.  The 
specimen  was  instrumented  with  one  25  mm  extensometer  and  two  piezoelectric  acoustic  emission 
sensors.  The  specimen  was  loaded  and  unloaded  six  times  (at  300  MPa,  400  MPa,  620  MPa,  710 
MPa,  845  MPa  and  failure)  such  that  each  loading  was  higher  in  magnitude  than  the  previous 
loadings.  Non-destructive  imaging  was  used  to  evaluate  the  damage  accumulation  after  each 
loading.  The  final  loading  caused  complete  failure  of  the  specimen.  A  low  displacement  rate  of 
0.127  mm/min  was  used  to  better  capture  acoustic  emission  events  as  a  function  of  time. 
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Tension  tests  with  fatigue  loading  were  conducted  on  an  MTS  universal  testing  machine 
according  to  ASTM  D3479  [4].  The  tests  were  conducted  on  three  quasi-isotropic  specimens  in 
order  to  thoroughly  characterize  the  quantity  and  location  of  damage  progression  as  a  function  of 
the  number  of  loading  cycles.  The  tests  were  conducted  at  a  loading  frequency  of  5  Hz,  an  r-ratio 
of  0.1,  and  a  maximum  stress  amplitude  of  143  MPa.  143  MPa  is  18%  of  the  maximum  strength 
(795  MPa)  of  an  independently  tested  quasi-isotropic  specimen. 

3.1.3  Acoustic  Emission 

Acoustic  emission  (AE)  testing  was  used  to  detect  failure  events  within  the  composite 
material.  In-situ  AE  activity  was  recorded  on  a  Micro-II  Digital  AE  System  produced  by  Physical 
Acoustics  Corporation.  When  a  material  experiences  local  failure,  it  releases  strain  energy  which 
produces  a  stress  wave  in  the  specimen.  The  AE  system  detects  this  acoustic  energy  and  records  it 
as  a  hit.  Prior  to  testing,  an  AE  calibration  study  was  performed  to  define  the  appropriate  signal 
conditioning  parameters.  It  was  found  that  an  amplitude  threshold  of  48  dB  enabled  the  detection 
of  all  valid  material  failure  events  without  recording  ambient  noise.  As  recommended  by  the 
equipment  manufacturer,  the  AE  timing  parameters  used  for  this  study  were  Peak  Definition  Time 
=  400  jus  ,  Hit  Definition  Time  =  800  jus  ,  Hit  Lockout  Time  =  200  jus  ,  and  Maximum  Duration 
=  100  ms. 

3.1.4  X-ray  Radiography 

The  monotonically  loaded  quasi-isotropic  specimen  was  loaded  to  ultimate  failure  in 
increments.  Once  each  load  level  was  achieved,  the  specimen  was  unloaded  and  examined  using  a 
160  kV  Philips  X-ray  system  (0.4  mm  focal  spot)  and  General  Electric  CR  Tower  Computed 
Radiography  system  using  IPS  imaging  plates  and  50  micron  sampling.  The  imaging  parameters 
were  26  kV,  3  inA,  and  30  s,  with  a  source-to-detector  distance  of  48  inches.  Prior  to  X-ray 
examination,  the  edges  of  the  specimen  were  exposed  to  zinc  iodide,  an  opaque  penetrant,  which 
was  absorbed  into  all  cracks  and  voids  adjacent  to  the  specimen  edge.  The  optimum  view  of 
damage  was  achieved  using  the  General  Electric  Rhythm  image  processing  software  where  a  level 
III  contrast  enhancement  filter  was  applied  to  each  X-ray  image  using  noise  reduction  and  latitude 
correction. 

X-ray  radiography  was  also  used  to  evaluate  the  state  of  damage  accumulation  for  the 
fatigue  loadings  at  1500,  25000,  50000,  and  100000  loading  cycles.  Once  the  chosen  number  of 
loading  cycles  was  reached,  the  specimen  was  removed  from  the  tensile  testing  machine  and 
examined  using  the  Philips  X-ray  system  with  the  previously  described  procedures  and  settings. 

3.1.5  X-ray  Computed  Tomography 

For  the  monotonically  loaded  specimens,  once  the  planar  X-ray  showed  a  significant 
amount  of  damage,  the  specimen  was  also  examined  using  an  X-Tek  HMX160  CT  system.  The 
main  components  included  an  X-ray  source,  a  rotation  stage  on  which  the  sample  was  fixed,  and  an 
X-ray  detector.  The  maximum  resolution,  at  highest  magnification,  was  approximately  5  /an .  A 
Molybdenum  target  was  used.  The  source  voltage  and  the  source  current  were  90  kV  and  90  juA , 
respectively.  The  specimen  was  clamped  vertically  approximately  33  cm  from  the  X-ray  source. 
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The  sample  translates  and  rotates  over  360°  with  a  step  size  of  0.5° .  Averages  of  eight  projection 
images  (1024  x  1024  pixels)  were  collected  at  each  position.  The  raw  image  data  was 
reconstructed  using  CT  Pro  software.  A  three  dimensional  structure  of  the  damaged  specimen  was 
visualized  in  order  to  evaluate  the  damage  through  the  thickness  of  the  specimen  using  3D  surface 
rendering  techniques. 

Acoustic  emission  testing  (AE),  X-ray  radiography,  and  X-ray  computed  tomography 
(CT),  are  instrumental  in  characterizing  some  aspects  of  damage  progression  and  model 
validation.  AE  uses  piezoelectric  sensors  to  passively  detect  acoustic  signals  emitted  by  the 
material  during  damage  propagation  [48,  33].  The  most  advantageous  characteristic  of  AE  is  that 
the  sensors  detect  damage  during  testing  in  a  range  that  cannot  be  distinguished  by  typical 
instrumentation  such  as  load  cells,  strain  gauges,  and  displacement  transducers.  X-ray  radiography 
is  a  common  NDI  method,  in  which  a  two  dimensional  image  is  recorded  on  an  imaging  plate  as 
energy  is  passed  through  a  stationary  material  [55,  60].  In-plane  delaminations  can  be  easily 
detected  due  to  the  variation  of  X-ray  absorption  between  the  material  and  the  void.  The  difficulty 
with  X-ray  radiography  is  the  inability  to  characterize  damage  as  a  function  of  specimen  thickness. 
X-ray  computed  tomography  provides  an  ultra-high  resolution  three  dimensional  image  through 
the  thickness  of  a  material  [56].  As  X-ray  CT  equipment  has  become  more  readily  available,  this 
technique  is  being  used  for  nondestructive  evaluation  of  composites  [53].  The  primary  advantage 
of  X-ray  CT  for  composite  materials  is  that  delaminations,  transverse  matrix  cracks,  and  fiber 
fracture  can  all  be  adequately  characterized  [18,  58]. 

3.2  Symmetric  Reduced  Order  Computational  Homogenization 

Mathematical  homogenization  theory  provides  a  rigorous  mathematical  framework  for 
modeling  the  response  of  heterogeneous  materials.  The  mathematical  theory  was  formalized  in  the 
seminal  works  of  Babuska  [6],  Bensoussan  [10],  Sanchez-Palencia  [51]  and  Suquet  [57],  among 
others.  Since  the  development  of  the  computational  framework  for  the  mathematical 
homogenization  theory  by  Guedes  and  Kikuchi  [32],  numerous  models  based  on  the 
computational  homogenization  method  (CHM)  have  been  proposed  to  predict  the  elastic  and 
inelastic  response  of  heterogeneous  materials  including  material  failure. 

The  distinct  feature  of  the  computational  homogenization  method  in  modeling  the  response 
of  heterogeneous  materials  is  in  the  evaluation  of  the  constitutive  response  at  a  material  point  of  a 
macroscopic  (homogenized)  medium.  In  CHM,  the  constitutive  response  of  the  equivalent 
homogeneous  medium  is  evaluated  by  solving  a  microscale  boundary  value  problem  defined  on  a 
representative  volume  element  (RVE)  of  the  heterogeneous  microstructure.  This  approach 
decouples  the  effect  of  the  microstructural  topology  from  the  material  behavior  of  the 
microconstituents,  as  well  as  the  conditions  along  the  microconstituent  interfaces.  CHM  simplifies 
the  constitutive  modeling  process  since  the  response  of  the  microconstituents  tend  to  be  simpler  to 
model,  compared  to  phenomenological  modeling  of  the  combined  microstructure-material 
behavior  effects.  In  the  case  of  modeling  the  failure  of  heterogeneous  materials,  a  number  of 
outstanding  computational  issues  remain,  including  selection  of  the  boundary  conditions  for  the 
RVE  problem  in  the  presence  of  defects  [62,  21],  evolution  of  the  RVE  domain  upon  defect 
formation,  size  scale  effects  [27],  and  spurious  mesh  dependency  [9],  among  others. 

One  additional  major  challenge  associated  with  the  computational  homogenization  method 
is  the  computational  cost  associated  with  solving  nonlinear  RVE  problems  to  evaluate  the 
constitutive  response  of  the  macroscopic  problem.  This  problem  is  alleviated  by  one  or  a 
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combination  of  two  approaches.  The  first  is  the  brute-force  parallelization  of  the  multiscale 
problem,  in  which,  the  RVE  problem  evaluations  are  distributed  to  a  large  number  of  compute 
nodes  and  evaluated  in  parallel  [20].  The  second  approach  is  reduced-order  evaluation  of  the  RVE 
problem.  Fast  Fourier  transfonn  [39],  proper  orthogonal  decomposition  [63],  spectral  method  [1], 
boundary  element  method  [30],  network  approximation  method  [11],  and  transfonnation  field 
analysis  (TFA)  [19,  7],  and  other  TFA-based  computational  methods  [24,  15,  38]  have  been 
effective  in  evaluating  the  inelastic  response  at  the  RVE  level  in  a  computationally  efficient 
manner.  In  a  recent  study,  eigendeformation-based  homogenization  method  (EHM)  was  proposed 
[45]  to  efficiently  evaluate  the  RVE  level  response  using  a  meso-mechanical  model.  This  method 
is  derived  based  on  a  generalization  of  the  transfonnation  field  analysis.  By  this  approach,  it  is 
possible  to  account  for  the  interfacial  debonding  effects,  in  addition  to  nonlinear  and  failure 
processes  within  the  constituent  materials  of  the  heterogeneous  microstructure. 

This  section  provides  a  model  reduction  methodology  for  efficient  evaluation  of  the 
microscale  boundary  value  problems  of  the  computational  homogenization  method.  The  presented 
approach  addresses  three  of  the  main  shortcomings  of  the  TFA-based  model  reduction  methods 
with  the  following  novel  contributions: 

1 .  A  new  methodology  for  the  determination  of  the  order  of  the  reduced  model  is  presented: 
The  accuracy  and  efficiency  of  the  reduced  models  clearly  depend  on  their  order  and  ability 
to  represent  the  failure  modes  within  the  microstructure.  A  reduced-order  model 
development  strategy  is  devised  to  identify  the  model  order  and  the  associated  coarse 
graining  at  the  microscale  for  accurate  and  efficient  representation  of  the  failure  modes. 

2.  The  proposed  reduced  order  model  leads  to  a  symmetric  formulation:  In  the  presence  of 
interfacial  debonding,  previous  eigendeformation-based  homogenization  formulations  lack 
symmetry,  which  increases  computational  cost. 

3.  The  proposed  formulation  eliminates  the  spurious  residual  stress  effect  upon  failure  due  to 
the  coarse  representation  of  the  inelastic  fields.  Some  of  the  transformation  field  analysis 
based  reduced  order  models  (e.g.,  [24,  45])  lead  to  spurious  residual  stress  fields  upon 
failure  in  the  microscale.  The  spurious  residual  stress  fields  pollute  the  macroscale  problem 
by  affecting  local  stress  redistributions. 

The  proposed  reduced  order  methodology  is  implemented  to  model  the  failure  response  of  brittle 
composite  systems,  in  which  the  failure  is  characterized  by  matrix  microcracking,  delamination 
and  debonding. 

The  remainder  of  this  section  is  organized  as  follows:  The  statement  of  the  multiscale 
problem  and  the  associated  macroscopic  and  microscopic  boundary  value  problems  are  presented 
in  Section  3.2.1.  In  Section  3.2.2,  formulation  of  the  symmetric  reduced  order  model  for  the 
microscale  problem  is  provided.  The  computational  algorithms  employed  to  evaluate  the  nonlinear 
reduced  order  model  are  discussed  in  Section  3.3.3.  Section  3.3.4  provides  small  scale  and  large 
scale  numerical  verification  examples  conducted  on  a  fiber  reinforced  matrix  composite. 
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3.2.1  Problem  Setting 


In  this  section,  we  present  a  summary  of  the  microscopic  and  macroscopic  boundary  value 
problems  associated  with  the  two-scale  asymptotic  homogenization  method  for  failure  response  of 
a  heterogeneous  body.  The  details  of  two-scale  asymptotic  homogenization  in  the  presence  of 
inelastic  effects  are  reported  in  the  literature  (see  e.g.,  Refs.  [61]). 

The  problem  setting  and  the  multiscale  heterogeneous  body  is  illustrated  in  Fig.  1.  The 
heterogeneous  domain,  denoted  by  Q ,  is  parameterized  by  the  macroscopic  coordinate  vector,  x. 
Q  is  composed  of  the  repetition  of  a  small  representative  volume  element,  0  ,  which  is 
parameterized  by  the  microscopic  coordinate  vector,  y .  The  size  scale  ratio,  ^  ,  between  the 
characteristic  lengths  of  the  representative  volume  element,  0  ,  and  the  macroscopic  body,  Q  is 
assumed  to  be  very  small,  such  that  a  first  order  asymptotic  decomposition  of  the  displacement 
field  is  sufficient  to  accurately  capture  the  response  of  the  material.  The  response  fields  are 
assumed  to  be  periodic  about  the  representative  volume  element.  The  periodicity  condition  states 
that  the  value  of  the  response  fields  are  the  same  at  the  opposing  faces  of  a  parallelepiped  RVE 
domain. 


Figure  1.  Macro-  and  Microscopic  Scales 


The  following  notation  is  employed  throughout  the  section,  unless  otherwise  noted: 
Subscript  roman  indices  denote  1,  2,  or  3.  Einstein  summation  convention  is  adopted  for  repeated 
indices.  Subscripts  xt  and  yi  following  a  comma  denote  differentiation  with  respect  to  the 

macroscopic  and  microscopic  coordinate  vectors,  respectively.  Differentiation  within  parentheses 
denotes  symmetric  differentiation  with  respect  to  the  indices.  Bold  characters  denote  tensor 
notation.  Macaulay  brackets  denote  averaging  over  the  RVE: 
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where,  |  0 1  is  the  volume  of  the  RVE. 

The  displacement  field  of  the  heterogeneous  body  is  expressed  using  a  two-scale 
asymptotic  expansion: 

«,■  (x,  y ,  t)  =  ut  (\,t)+ £u]  (x,  y ,  t)  (2) 

in  which,  u  is  the  macroscopic  displacement  field,  and;  u1  is  the  variation  of  the  displacement 
field  within  the  RVE. 

Microscale  Problem 


In  the  presence  of  failure  processes,  u1  is  described  by  the  microscopic  equilibrium 
equation  defined  over  the  RVE  (i.e.,  ye0) 

\ju  (y )  fa  M+ u{k,y,)(x’  y,*)~  Ma  (x>  y  >  *)] 


k 


=  0 


(3) 


in  which,  L  is  the  fourth  order  tensor  of  elastic  moduli,  taken  to  be  symmetric  and  strongly 
elliptic,  s  =  V' u  the  macroscopic  strain  tensor;  V*(-)  =  (-)(j  x  )  denotes  the  symmetric  gradient 

operation  with  respect  to  macroscopic  coordinates;  and  //  damage  induced  inelastic  strains.  In 
this  work,  the  damage  induced  inelastic  strains  are  modeled  using  a  scalar  continuous  damage 
mechanics  model: 

Mu  (x,  y,  t)  =  coph  (x,  y,  t)sy  (x,  y,  t)  (4) 


in  which,  coph  e  [0,1 )  is  a  history  dependent  variable,  which  represents  damage  within  the 
microconstituents,  and  e  is  the  strain  tensor.  Using  the  scaling  relations  provided  by  the 


asymptotic  decompositions  with  multiple  spatial  scales: 

£ij (x, y, t)  =  Sj{x,t)+  ul-.  t  j (x, y, t)  (5) 

Along  the  microconstituent  interfaces,  debonding  is  considered  based  on 
traction-separation  laws  given  as  ( y  e  5” ) 

^ (x, y , i) - [l - coint (x, y , t )]k jV (y )SN (x, y , t) < 0 ;  AA'(x,y,t)>0  (6) 

{fv(x,y,t)-[l-n>,,,(x,y,t)]kw(y)AiV(x,y,t)}AjV(x,y,t)  =  0  (7) 

tr(x,  y,  t)  =  [l  ■ -  ®,„(x,  y,  t)]kT  (y)c7r  (x,  y,  t)  (8) 


in  which,  0)int  e  [0,1)  is  a  history  dependent  variable,  which  represents  damage  along  the 

interface;  t A  and  8  '  are  the  components  of  the  traction  and  displacement  jump  normal  to  the 
interface,  respectively;  k  '  (y)  and  k 1  (y )  the  initial  interface  stiffness  in  the  normal  and 

tangential  directions,  respectively,  and;  tT ,  8  r  the  tangential  components  of  the  traction  and 
displacement  jump  along  the  interface,  respectively.  The  traction  and  displacement  jump 
components  are  expressed  in  tenns  of  the  local  coordinate  system  formed  by  the  normal  and 
tangential  directions  at  the  interface  point. 

The  microscale  problem,  which  is  a  nonlinear  boundary  value  problem  is  solved  to 
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evaluate  the  microscale  displacement  field  u1  by  imposing  periodic  boundary  conditions  along 
the  exterior  boundaries  of  the  RVE  while  restricting  the  rigid  body  motion.  The  microscale 
boundary  value  problem  is  quasi-static  as  indicated  by  the  lack  of  inertial  terms  in  the  governing 
equations.  The  present  fonnulation  is  limited  to  the  cases  for  which  the  characteristic  size  of  the 
RVE  is  small  compared  to  the  length  of  the  deformation  and  stress  waves. 

Macroscale  Problem 


The  macroscopic  displacement  field  is  described  by  the  macroscopic  momentum  balance 
equation  defined  over  Q  : 

Vy,x  .{x,t)+bi{x,t)  =  p(x,  t%  (x,  t)  (9) 

J 

in  which,  double  dot  over  a  field  denotes  twice  differentiation  in  time;  a  denotes  the 
macroscopic  stress  tensor,  evaluated  by  volume  averaging  of  the  stresses  over  the  domain  of  the 
RVE 

ov(x,t)=(cr..)  (10) 

The  stress  field  is  expressed  as: 

Vy  (x,  y,t)=  Lm  (y  )[skt  (x,  t)  +  u\kyi)  (x,  y,  t)  -  nkl  (x,  y ,  f  )J  (11) 

b  and  p  denote  the  RVE-average  body  force/unit  volume  and  the  RVE-average  density, 
respectively: 

bi(x,t)=(bi);  p  =  (p)  (12) 

The  boundary  and  initial  conditions  of  the  macroscale  initial-boundary  value  problem  are 
defined  as 

ui(x,t)  =  n,(x);  xeQ;  t  =  0  (13) 

uj(x,t)  =  v,(x);  xeQ;  t  =  0  (14) 

ui(x,t)=ui(x,t);  x  e  TI( ;  t  e  [0, t0 ]  (15) 

vij(x,t)nj=ti(x,t\  xer,;  le[0,io]  (16) 

in  which,  u ,  u  are  prescribed  initial  and  boundary  displacements,  respectively;  v  prescribed 
initial  velocity,  and;  t  prescribed  boundary  traction.  The  prescribed  initial  and  boundary 
conditions  are  assumed  to  be  constant  with  respect  to  the  microscopic  coordinate  vector  y . 

3.2.2  Reduced  Order  Modeling  of  the  Microscale  Problem 

The  macroscale  problem  defined  in  Section  3.2.1  is  coupled  with  the  microscale  problem 
defined  in  the  same  section  through  the  macroscopic  constitutive  relationship  (Eqs.  10  and  1 1). 

The  evaluation  of  the  macroscopic  stress  at  each  macroscopic  material  point  requires  the  solution 
of  the  microscopic  RVE  problem  associated  with  that  material  point.  When  the  finite  element 
method  is  employed  to  evaluate  the  macroscale  problem,  a  nonlinear  microscale  problem  must  be 
evaluated  to  update  the  stress  at  each  integration  point  for  each  increment  and  iteration  of  every 
time  step  of  the  loading  history.  This  is  a  tremendous  computational  burden.  In  this  section,  a 
novel  reduced  order  model  is  derived  to  efficiently  compute  the  microscopic  response.  To  this 
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extent,  the  microscale  displacement  field  is  decomposed  into  linear  and  damage  induced 
components: 

u)  (x,  y  ,t)  =  Hikl  (y  )skl  (x,  t )  ■ +  u,  (x,  y  ,t)  (17) 

in  which,  H  is  the  third  order  elastic  influence  function  obtained  by  substituting  Eq.  17  into  Eq.  3 
and  solving  the  microscale  problem  in  the  absence  of  all  inelastic  processes  (i.e.,  coph  =  0)int  =  0). 

u  is  the  displacement  field  induced  by  the  damage  processes  within  the  microconstituents  and  the 
interface: 

«/ (x, y , t)  =  [h’t,  (y , y K (x, y , t)dy  +  (y , y fe, (x, y,/>/y  (18) 

in  which  hp/'  and  h"“  are  the  phase  damage  and  interface  damage  induced  influence  functions. 
h/,/!  and  h"!'  are  the  particular  solutions  to  the  RVE  problems  obtained  by  substituting  Eq.  17 
into  Eq.  3  and  solving  the  microscale  problem  in  the  presence  of  phase  damage  (i.e.,  // )  and 
interface  damage  (i.e.,  8),  respectively.  The  governing  equations  and  the  discrete  approximations 
of  the  elastic  and  damage  induced  influence  functions  are  provided  in  Ref.  [45]  and  will  not  be 
discussed  herein.  In  this  section,  we  concentrate  on  the  new  model  reduction  methodology  based 
on  the  microscopic  displacement  field  decomposition  provided  in  Eqs.  17  and  18. 

Substituting  Eq.  17  into  Eq.  3,  premultiplying  the  resulting  equation  with  hph ,  and 
integrating  over  the  domain  of  the  RVE  yields: 

(v  y  )fcy™  (y)[Amnki  (y  )%  (x,?)+  ?m„  (x,  y,  t)-/Jmn  (x,  y,  ?)]}^  dy  =  0  (19) 

in  which,  a  =  V'u  ;  A  =  I  +  G  is  the  fourth  order  elastic  strain  concentration  tensor;  I  the 

fourth  order  identity  tensor,  and;  G  =  V*H  ,  the  elastic  polarization  tensor.  The  use  of  Eq.  19 

secures  a  symmetric  formulation  as  subsequently  derived.  This  is  in  contrast  with  the  previous 
eigendefonnation-based  reduced  order  models,  which  are  non-symmetric  [45].  Integrating  by 
parts,  applying  divergence  theorem  and  employing  the  perodicity  of  the  response  fields  over  the 
domain  of  the  RVE  yields: 

J/^(y»y)^(y)[4»^(y)^(^0+^»»(x,y,0-//m»(x>y^)]^  =  0  (2°) 

where,  gph  =  V(h,,/'  is  the  fourth  order  phase  damage  polarization  tensor. 

A  second  set  of  equilibrium  equations  are  obtained  by  premultiplying  the  microscale 
equilibrium  equation  (Eq.  3)  with  \\mt ,  and  following  a  similar  procedure  as  described  above: 

\&sZ  (21) 

in  which,  gint  =  Vsyti'“  is  the  third  order  interface  damage  polarization  tensor.  Substituting  Eq.  4 


into  Eqs.  20  and  21  yields: 

JJi  -  a#  (v  y,  OkS  (y>y  (y  )• 

•  Umnk,  (y  K/  (v  t)  +  smn  (x,  y,  t)]dy  =  0 

l^ph  (y,y  )Lijmn  (y)- 

•  Umnki  (y)%  (v  t )  +  smn  (x,  y,  t)\ly  =  -tp  (x,y,  t) 


(22) 

(23) 


We  introduce  the  following  discretizations  for  damage  fields  coph  and  cojnl ,  and  the 
damage  induced  fields,  //  and  8  using  mesomechanical  shape  functions 
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®ph >  Mu  Xx»  y>  *) = Z#J*(y  )H^  /4r)  )(x,  t) 


fain,  A  J(x,y,0= EAr^)(y)H»,’^(/?)J(x’r)  (25) 

p=\ 

in  which,  y  —  1 ,2, . . . ,  n  and  /?  =  1 ,2, . . . ,  in  ;  «  and  m  denote  the  level  of  discretization  within 
the  phases  and  along  the  interface,  respectively.  8  denotes  the  displacement  jump  vector  in  the 
local  coordinate  system  (i.e.,  8 (/?)  =  ST^'\  ).  The  phase,  N^h\  and  interface,  N\^ ,  shape 

functions  have  compact  support  within  subdomains  of  the  phases  and  the  interface 

7VW(y)=O?/y^0(/);  0<r)  <=  0  (26) 

N^(y)  =  0ify£SW-  S(P)aS  (27) 

Employing  Eqs.  24  and  25,  s  is  expressed  in  terms  of  the  damage  induced  strain  and 
displacement  jump  coefficients: 

s,(w  )= m 

r  P 

in  which,  the  coefficient  tensors  P  and  R  are: 

Pm  (y ) = \e(r)si  (y»  y  W’  (y  Vy  (29) 

^f(y) = \s(P)sTq  (y>  yK?(y)e?  (y  Vy  (30) 

where  eq  denotes  the  transformation  vector  between  the  global  and  local  coordinate  systems 
along  the  interface. 

Substituting  Eqs.  24  and  25  into  Eq.  22,  premultiplying  the  resulting  equation  with 
N$(y),  and  integrating  over  the  domain  of  the  RVE  yields  (rj  —  1 ,2, . . . ,  n  ): 

it*-  ^t^X’OLA^t^yft^y^iy)  (3d 

A— 1  ^ 

[Anna  (y  (x,  t) + smn  (x,  y ,  tyfa } = o 

Similarly,  substituting  Eqs.  24  and  25  into  Eq.  23,  premultiplying  the  resulting  equation  with 
N^t\ y)  and  integrating  over  the  domain  of  the  RVE  yields  ( a  —  1 ,2, . . . ,  m  ): 

Z}1  -  ®S)(x,0]fe(A^S)(yft)(yXy™(y)  (32) 

[Amnk,  (y%  (x,  t) + y ,  t)}iy }  =  -if]  (x,  t) 

Substituting  Eq.  28  into  Eqs.  31  and  32  results  in 


m 

=o 
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d{& )£u  (x^  *) + vi-f j  (x^) 


r= i 


m 


(34) 


■--‘PM 


where, 


in  which,  C(,/Vi, 
coefficient,  t'a> 


D 

is: 


(aA) 


=  Jq(a®  (y)Lqrmn(y)\nkh)NtAy)dy 

(35) 

=  \&{^\y)hmn(y)Annkh  K^y^y 

(36) 

F^r) = [(ASlyfvJyKfyW’iyfe 

(37) 

Jtp) = l0(A)^(yK,™«(y)^!S(y  K?(y  Vy 

(38) 

m(mp) = £(A)^)(y)4-m«(y®!(y  Wt}(y)^y 

(39) 

('/A;/| ,  Ji:,/  V'1  and  M!"v;|  are  coefficient  tensors.  The  interface  traction 

(40) 

The  relationship  between  the  interface  traction  and  the  displacement  jump  is  nonlinear.  It 
is,  therefore,  not  possible  to  derive  explicit  expressions  for  the  relationship  between  the  traction 
and  displacement  jump  coefficients.  In  this  section,  the  relationship  between  the  pointwise 
tractions  and  displacement  jumps  are  adopted  to  represent  the  relationship  between  the  traction 
and  displacement  jump  coefficients.  This  approach  has  also  been  employed  in  a  number  of 
previous  investigations  (e.g.,  [38]).  The  unilateral  contact  and  adhesion  conditions  are  expressed 
as 

^(o)(x,r)>0  (41) 

^(a)(x,0-[l-®£)(x.0Ka)^JV(“)(x.0^0  (42) 

0  (43) 

The  tangential  adhesion  condition  is  also  written  in  a  similar  form  as 

(44) 

Similar  to  the  interface  traction-separation  conditions,  the  nonlinear  evolution  of  the  phase 

,i  '/  l  are  expressed  in  terms  of  the  field  coefficients: 

(45) 

(46) 

in  which,  q  and  qmt  are  state  variables  defining  the  evolution  of  the  phase  and  interface  damage 
variables,  respectively,  and; 

[•I  ;/’(x,  t)  =  \e(y)N{;2  (y)[-](x,  y,  t]dy-  [•]'*’  (x,  t)  =  {^^(yHx,  y,  t)dy  (47) 

Equilibrium  equations  (Eqs.  33  and  34),  in  addition  to  the  interface  conditions  provided  by 
Eqs  41-44,  and  the  evolution  equations  for  the  phase  and  interface  damage  coefficients  form  the 


Ph  and  col,  > 


and  interface  damage  coefficients,  aiy) 

,l(o 

k,s„€') 


Mph  -  °>Ph 


co  —  CO  , , 

int  int  \ 


->  coya  =  co[at 

int  int 


(/(«)  x(a)  nint(a)' 
X  i  ’ui  ’Oa 
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reduced  order  model.  The  reduced  order  model  is  solved  to  obtain  the  unknown  coefficients  py) 
and  S(a). 


The  macroscopic  stress  tensor  is  expressed  in  terms  of  p(r) ,  Sla> ,  and  the  macroscopic 
strain  tensor  by  substituting  Eq.  17  and  18  into  Eq.  1 1  and  using  Eqs.  24  and  25  to  obtain 


where,  the  coefficient  tensors  L 


TtA)  p(Ar) 


r 

r(A/9) 


(48) 


and  R  7  are  expressed  as 


4*  =  7^7  l^K'(y)L,Jy)A,Jy)dy  (49) 

P,l if1  =  ^  1; » ,^7  fj'  y  )?;i  (y  k/y  (50) 

=7^7jewA'Sl(y)is.,(yR(£(yhy  (si) 


3.2.3  Computational  Aspects 

The  implementation  of  the  proposed  reduced-order  multiscale  model  is  conducted  in  two 
stages.  The  preprocessing  stage  consists  of  detennining  the  model  order,  partitioning  of  the  RVE 
based  upon  the  model  order,  and  computing  the  coefficient  tensors  associated  with  the  reduced 
order  model.  The  macroscopic  analysis  stage  consists  of  evaluating  the  macroscale  problem 
described  in  Section  3.2.1  using  a  numerical  method.  In  this  study,  the  macroscopic  analyses  are 
conducted  using  the  finite  element  method.  The  commercially  available  finite  element  analysis 
program,  Abaqus,  is  employed.  Reduced  order  multiscale  models  and  direct  numerical  simulations 
for  the  verification  of  the  proposed  approach  are  conducted  using  the  user  supplied  subroutine 
capabilities.  The  remainder  of  this  section  discusses  a  novel  strategy  for  selection  of  the  model 
order  and  partitioning  of  the  RVE  domain,  as  well  as  the  numerical  evaluation  of  the  reduced  order 
model. 

Reduced-Order  Model  Development  Strategy 

The  proper  partitioning  of  the  RVE  domain  is  critical  to  the  efficiency  and  the  accuracy  of 
the  proposed  reduced  order  modeling  approach.  The  partitioning  of  the  RVE  consists  of  the 
selection  of  the  number  of  phase  (  n  )  and  interface  ( m  )  partitions,  as  well  as  the  domain  of  each 
phase  ( 0(r) )  and  interface  ( Slfl> )  partition.  Theoretically,  as  the  number  of  partitions,  n  and  m  , 
increase,  the  accuracy  of  the  reduced  model  increases  at  the  expense  of  additional  computational 
cost.  The  accuracy  of  the  reduced  order  model  is  also  strongly  affected  by  the  selection  of  the 
partition  domains,  0(/l  and  Si/]>  for  a  given  number  of  phase  and  interface  partitions.  Previous 
and  current  investigations  found  significant  variability  in  the  performance  of  a  reduced  order 
model  based  on  the  partitioning  strategy.  It  is  possible  to  adaptively  select  the  order  of  the 
reduced-order  model  based  on  a-priori  error  measures  associated  with  the  state  of  the  failure 
processes  during  the  macroscopic  simulations.  Such  a  dynamic  partitioning  strategy  was  proposed 
in  Ref.  [45],  The  dynamic  strategy  resembles  h  -version  adaptive  finite  element  modeling  with 
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goal  oriented  mesh  adaptivity  [42],  or  multilevel-multiscale  modeling,  in  which  the  heterogeneity 
within  the  process  zones  are  adaptively  resolved  [28].  The  dynamic  partitioning  strategy,  while 
rigorous,  comes  with  a  significant  increase  in  computational  cost.  The  additional  computational 
cost  is  due  to  error  assessment  and  recomputation  of  the  coefficient  tensors  during  the  macroscopic 
analysis. 


ttttttt  tt 
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Figure  2.  The  Partitioning  and  Model  Reduction  Strategy 


In  this  section,  a  novel  static  partitioning  strategy  is  presented,  in  which,  the  RVE  domain 
partitions  and  the  model  order  is  identified  prior  to  the  macroscopic  analysis.  The  present  approach 
provides  a  model  selection  strategy  capable  of  accounting  for  the  failure  modes  within  the 
microstructure  using  a  small  number  of  domain  partitions.  The  model  selection  strategy  consists  of 
identification  of  the  failure  paths  within  the  microstructure  when  subjected  to  a  number  of  loading 
modes,  and  partitioning  the  domain  of  the  RVE  as  well  as  the  interfaces  by  selecting  each  failure 
path  as  a  partition.  The  failure  paths  within  the  microstructure  are  identified  by  conducting  detailed 
RVE-level  simulations.  The  RVE  is  subjected  to  uniform  macroscopic  strain  modes  (e.g.,  uniaxial 
tensile  or  compression  and  shear).  Figure  2  illustrates  the  identification  of  the  failure  paths  in  a  2-D 
particle  reinforced  matrix  under  uniform  macroscopic  axial  and  shear  strains.  The  failure  path  due 
to  each  loading  condition  is  marked  as  an  individual  partition  as  shown  in  Fig.  2.  RVE-level 
simulations  were  conducted  by  applying  Dirichlet  conditions  along  the  boundaries  when 
determining  the  failure  paths.  This,  along  with  the  unstructured  finite  element  mesh  leads  to 
unsymmetric  failure  paths  despite  symmetry  in  the  RVE  geometry.  Periodic  boundary  conditions 
are  maintained  along  the  RVE  boundaries  in  the  multiscale  model. 

The  failure  paths  associated  with  different  loading  modes  intersect  each  other  as 
demonstrated  in  Fig.  2.  Hence,  the  phase  partitions  are  allowed  to  overlap.  The  phase  shape 
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functions,  are  selected  to  accommodate  such  an  overlap.  Let  the  domain  of  the  RVE  be  partitioned 
into  n  possibly  intersecting  subdomains,  denoted  by  0(,y) ,  y  =  1,2,...,  n  .  The  union  of  the 
subdomains  spans  the  domain  of  the  RVE: 

0  ss  |j0(r)  (52) 

7=1 

The  intersection  between  two  partitions  are  denoted  as  0(;/7?'1  =  0(;/)  n  0(,,) .  A  material  point 
within  the  RVE  may  lie  in  all  n  partitions  or  less.  The  intersections  between  multiple  partitions 

are  defined  by  repetitive  Greek  superscripts:  0(y'?v"')  =  0(;/)  n0(,,)  n  0(r) _ The  shape  functions 

for  the  reduced  order  model,  N{^ ,  are  chosen  as: 


in  which,  i  —  ,  and; 


1 

i 

0 


if  ye  0|7) 
elsewhere 


(53) 


&[r)  =0w\|j0 


(in) 


n= i 


@;rt=U 

7J=\ 


0(w)\|j0 


{yr/v) 


v=l 


(54) 

(55) 


The  expressions  for  ©3 . . .  0;’  are  derived  analogously.  The  shape  functions  defined  in  Eq.  53 
allow  the  possibility  of  intersecting  shape  functions,  and  satisfy  the  partition  of  unity  property: 


2X*(y)=1;  ye®  (56) 

r= 1 

The  interface  shape  functions  are  continuous  across  the  partitions  to  satisfy  the  continuity 
of  tractions  and  displacement  jumps  across  the  interface  partitions.  Consider  the  partitioning  of  the 
interface  S  into  m  overlapping  subdomains  5,").  The  interface  shape  function,  Vi(j() ,  is  a 
linear  combinations  of  standard  finite  element  shape  functions  corresponding  to  the  nodes  along 
the  interface  partition,  S{a)  [45] 

vifiyE  TnM  jes  (57) 


in  which,  Na  is  the  standard  finite  element  shape  function  associated  with  the  microscopic  finite 
element  mesh  node  a  . 


Numerical  Evaluation  of  the  Reduced-Order  Model 


The  evaluation  of  the  reduced  order  model  for  the  microscale  problem  constitutes  the 
macroscopic  stress  update  at  a  macroscopic  point.  The  reduced  order  model  is  evaluated  using  the 
active  set  strategy  to  account  for  the  contact  conditions  at  the  interfaces. 

Given:  At  a  macroscopic  material  point  x  and  at  time  t ,  the  equilibrium  state  defined  by 
the  macroscopic  strain  tensor  ts  ;  the  inelastic  strain  and  displacement  jump  coefficients,  tju(r) 

and  ,S(a) ,  respectively,  where  y  =  and  a  =  1 ,2, ..., m  ;  state  variables  ,q<;  )  and 


16 

Approved  for  public  release;  distribution  unlimited 


t(\int(a) ,  which  define  the  evolution  of  the  phase  and  interface  damage  state,  ,co(^  and  tcof?  , 

respectively;  as  well  as  the  change  in  the  macroscopic  strain  state,  As  (taking  an  assumed  strain 
approach  in  the  numerical  evaluation  of  the  macroscopic  boundary  value  problem). 

Compute:  The  current  values  (at  time:  t  +  At )  of  the  inelastic  strain  and  displacement 

coefficients,  p  f)  and  5(a) ,  respectively;  the  current  damage  state,  «Kh  and  state 


variables,  q(r)  and  q.n,(a°  and  the  macroscopic  stress,  a  . 

In  this  section,  we  will  employ  vector  notation  using  the  classical  index  contractions  (e.g., 

1  =  {Lij})  \Lyki  }).  A  left  subscript  t  indicates  the  value  of  the  function  at  time  t.  A  left 
superscript  denotes  iteration  count.  The  eigendeformation  vector  is  defined  as: 

d  =  {//(l),...,//(”U(1),..^(m)f  (58) 

The  active  set  is  defined  as  the  set  of  all  interface  partitions  in  which  normal  displacement  jump 
coefficients  are  zero: 


A  =  \a  \  a  e  {l 

We  define  the  discrete  system  of  nonlinear  equations, 

in  which, 

};  SN(a)  =  0) 

VF ,  based  on  reduced  order  model  as: 

(59) 

(60) 

A  =1 

K(A) 

(61) 

and, 

1 

> 

> 

a 

g(1a1)  ••• 

£j(lAm) 

K^  = 

j^(wA0  p(«A«) 

^j(lAl)  ^  ^j(lAw) 

^(«Al) 

H(1A1)  ••• 

^-^(«Am) 

jj(lAm) 

(62) 

jj(mAl) 

jj(mAm) 

K 1 A 1  are  symmetric  matrices  since  G(“A'/  I  = 

F(^)=Fr(,A,)  and 

_  jjr(/?Aa) 

^ trac 

is  defined  as 


Kt 


1 

•  •  O 

...  0 

•  •  o 

1 

O  •• 

•  o 

...  o 

0 

...  o 

•  o 

...  0 

k« 

...  0 

1 

o  •  • 

...  0 

0 

1 

Jl 

<-* 

k(a)  =  (l-«;(or) 


int  I 


kf>  0 


d«) 

■N 

0 

0 


0 

kf]  0 


d«) 

VT 

0 


M) 


(63) 


The  symmetry  of  K  leads  to  this  formulation  being  denoted  as  a  symmetric  formulation.  The 
solution  of  T  =  0  is  evaluated  by  symmetric  nonlinear  solvers  rather  than  unsymmetric  ones  as 
has  been  the  case  in  previous  fonnulations.  f  is  defined  as 


A=1 
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and, 

f(A)  =  {c(1A),...,C("A),D(1A),...,D(mA)}r£  (65) 

In  the  case  of  tensile  loading  throughout  the  interface,  the  active  set  is  empty  (A  =  0  ), 
y(d)=0  solves  the  reduced  order  model.  When  A  ^  0  a  reduced  system  of  equations  is 
defined: 

>r4(d)=KAdA+fA  (66) 

in  which,  dA  and  fA  is  constructed  by  removing  each  row  which  corresponds  to  for  each 
partition  a  in  A  ,  from  d  and  f ,  respectively.  KA  is  constructed  by  removing  each  row  and 
column,  which  corresponds  to  8^)  for  each  partition  a  in  A  ,  from  K  .  The  rows  removed 
from  Eq.  60  form: 

1,l(<l)=KA<lA+fS  («7) 

The  reduced  order  model  is  evaluated  by  ensuring  TA(d)=  0  and  vPA(d)>0  are  satisfied.  The 
latter  condition  is  necessary  to  ensure  negative  tractions  upon  compressive  loading  along  the 
interfaces.  The  computational  algorithm  to  evaluate  the  reduced  order  problem  based  on  the  active 
set  strategy  is  provided  in  Box  2.  The  algorithm  is  initiated  by  setting  the  working  set,  which 
approximates  the  active  set  at  ( t  +  At ),  to  the  active  set  at  time  t ,  as  well  as  setting  the 
eigendefonnation  vector  to  d  (Step  1).  Within  each  iteration  k  ,  the  test  eigendeformation 

vector,  k  d  is  evaluated  by  a  standard  nonlinear  root  finding  algorithm,  such  as  Newton-Raphson 
or  quasi-Newton  methods  (Step  2a).  In  this  section,  a  quasi-Newton  SRI  method  [41]  is  employed 
to  compute  the  roots  of  VF k  .  In  this  method,  a  symmetric-rank-one  matrix  is  added  to  an 

approximation  of  the  Jacobian  of  VP k  at  each  iteration  of  the  nonlinear  solver.  In  practice,  this 

update  has  been  shown  to  provide  very  good  approximations  of  the  Jacobian  resulting  in 
superlinear  convergence.  The  advantages  of  quasi-Newton  methods  are  that  they  do  not  require  an 
explicit  fonnula  for  the  Jacobian  and  that  the  update  can  be  performed  on  the  inverse  Jacobian 
alleviating  the  need  to  solve  a  linear  system  of  equations  at  each  iteration.  In  this  paper,  the 
algorithm  is  initialized  with  a  finite  difference  approximation  to  the  Jacobian.  If  the  computed 
normal  displacement  jump  coefficients  violate  the  impenetrability  condition  (Steps  2b-c),  the 
partition  with  the  most  severe  violation  is  added  to  the  working  set. 

The  active  set  algorithm  for  evaluation  of  the  reduced  order  model  with  unilateral  contact 
constraints: 

1 .  Initialize  the  algorithm  by  setting  the  initial  guess  for  the  eigendeformation  vector  and  for  the 

active  set: 

k  =  l;  °d  =t  d;  *W=  A 

in  which,  k  W  denotes  the  working  set.  The  working  set  is  an  approximation  to  the  active 
set  at  iteration  k . 

2.  Loop  over  the  iterations  k : 
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(a)  Compute  d  by  evaluating  vFi.w(d  J=  0  using  a  standard  nonlinear  root  finding 

algorithm. 

(b)  Loop  over  each  interface  partition,  a  ,  which  is  not  in  the  working  set 
(i.e.,  a<£k\N): 

i.  If  the  impenetrability  condition  at  partition  a  is  violated  (i.e.,  k  ’  <  0 ), 

compute  the  step  size  0  <  A(a)  <  1  for  each  interface  partition  violating  the 
impenetrability  condition  as: 

k-lc(a) 

A(a>  =  — - - 

k  £(<*)  _k~l  £(a) 

UN  UN 

(c)  If  the  step  size  is  reduced  (i.e.,  0 ): 

i.  Compute  ft  =  argmin{/t<o0 } 

ii.  Update  the  working  set:  k+]  W-  Wu  {/?} 

iii.  k  ^ —  k  +  l 

iv.  td  =  i')(td-*-,d)+*-,d 

v.  Return  to  the  beginning  of  the  iteration  loop 

(d)  Check  if  the  unilateral  conditions  are  violated  in  any  partition  within  the  working  set: 

(i.e.,  If  any  component  of  vPA_(d)<  0) 

i.  Compute  J3  =  argminj1!7,.  ^  ) 

ii.  Update  the  working  set:  k+]  W  =k  W  -  {/?} , 

iii.  k  ^ —  k  +  \ 

iv.  *d=*d 

v.  Return  to  the  beginning  of  the  iteration  loop 

(e)  Update  the  eigendeformation  vector:  d  =/  d 

(f)  Update  the  active  set:  A  =k  W 

(g)  Exit  the  algorithm 
End  iteration  loop 

When  the  computed  normal  displacement  jump  coefficients  do  not  violate  the 
impenetrability  condition,  the  unilateral  contact  constraints  are  checked  within  the  active  set.  If  the 
unilateral  contact  conditions  are  violated  (i.e.,  if  the  computed  interface  traction  coefficients  are 
positive  at  partitions  within  the  working  set),  the  partition  with  the  most  severe  violation  (largest 
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positive  interface  traction  coefficient)  is  removed  from  the  active  set  (Step  2d).  When  the 
unilateral  contact  constraints  are  satisfied,  the  active  set,  the  eigendeformation  vector,  the 
associated  internal  state  variables,  and  damage  variables  are  updated  (Steps  2e-g). 

Two-Order  Reduced  Modeling 

Reduced-order  models  fail  to  accurately  capture  the  post-failure  response  of  the 
representative  volume  element.  The  failure  is  defined  as  the  loss  of  load  carrying  capacity  along  at 
least  one  loading  direction.  For  instance,  full  damage  within  any  one  of  the  failure  paths  along  with 
interface  debonding  in  the  RVE  illustrated  in  Fig.  2  cause  failure  along  the  associated  load 
direction.  The  reduced  order  models  exhibit  spurious  residual  stiffness  upon  failure,  which 
prohibits  proper  redistribution  of  the  stresses  at  the  macroscopic  scale.  While  increasing  the  model 
order  diminishes  the  spurious  residual  stiffness,  this  approach  increases  the  computational  cost. 

We  propose  a  two-order  modeling  scheme  to  eliminate  the  residual-stress  fields  upon 
failure  without  significantly  compromising  the  computational  efficiency.  In  this  approach,  the 
stresses  are  computed  based  on  a  high  order  model,  whereas  the  damage  coefficients  are  evaluated 
using  the  low-order  reduced-order  model  described  in  Section  3.2.2.  The  stress-update  procedure 
for  the  two-order  reduced  model  is  as  follows: 


1 .  Evaluate  the  eigendeformation  vector,  dlow ,  and  damage  coefficients,  coW ; 

7  =  1,2,...,  nlow  and  0)\" ) ;  a  =  1,2,...,  mlow  for  the  low-order  model  using  the 


numerical  procedure  described  previously  in  this  section.  nlow  and  mlow  are  the  orders 


for  the  low-order  model  selected  by  the  reduced-order  model  development  strategy 
described  in  Section  3.2.2. 

2.  Map  the  damage  coefficients  of  the  low-order  model  onto  the  high-order  model 
partitions.  The  mapping  of  the  damage  coefficients  onto  the  high-order  model  partitions 
is  trivial  when  the  high-order  model  is  constructed  by  hierarchical  subpartitioning  of  the 
low-order  model.  In  this  study,  each  finite  element  within  the  RVE  domain  constitutes  a 
partition  for  the  high-order  model. 

3.  Evaluate  the  eigendeformation  vector,  dhigh ,  for  the  high  order  model  by  solving  the 


linear  system: 


*  high 


—  J£-1  f 

^  high  A  high 


(68) 


4.  Compute  macroscopic  stress  (Eq.  48),  using  the  eigendeformation  vector  of  the 
high-order  model. 


3.2.4  Numerical  Verification 

The  capabilities  of  the  proposed  reduced  order  modeling  methodology  are  verified  against 
direct  finite  element  simulations.  The  verification  study  consists  of  (1)  analysis  of  an  RVE 
response  and  assessment  of  the  reduced  order  model  predictions,  and  (2)  a  three-point  bending 
problem  to  assess  the  capabilities  of  the  reduced  order  model  in  capturing  the  overall  failure 
response  of  macroscopic  structures. 
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RYE  Analysis 


The  multiscale  methodology  described  in  the  previous  sections  is  applied  to  develop  a 
meso-mechanical  model  for  a  2-D  composite  matrix  with  a  circular  inclusion.  Figure  2  illustrates 
the  geometry  of  the  microstructure.  The  evolution  of  damage  within  the  matrix  and  along  the 
interface  is  modeled  based  on  continuous  damage  mechanics  models  proposed  in  [45]  for  brittle 
composite  constituents.  The  reinforcement  is  assumed  to  behave  elastically  within  the  range  of 
applied  loads. 

The  capabilities  of  the  proposed  multiscale  model  in  capturing  the  failure  modes  for  a 
range  of  loading  conditions  are  verified  by  comparing  the  model  simulations  to  direct  numerical 
simulation  of  the  representative  volume  element.  The  finite  element  mesh  employed  in  these 
simulations  is  shown  in  Fig.  2.  The  characteristic  material  length  scale  associated  with  the  matrix 
constituent  is  assumed  to  be  1/8  of  the  RVE  size.  The  finite  element  mesh  of  the  RVE  is  designed 
to  have  an  average  size  of  1/8  of  the  RVE  length  scale  to  avoid  numerical  errors  associated  with 
mesh-sensitivity.  The  multiscale  model  is  developed  based  on  the  reduced  order  model 
development  strategy  described  in  Section  3.2.2.  The  reduced  order  model  is  developed  using  the 
biaxial  tension,  uniaxial  tension  and  shear  loading  modes.  The  partitioning  of  the  reduced  order 
model  is  shown  in  Fig.  2.  The  matrix  phase  and  the  interface  are  modeled  using  6  and  4  partitions, 
respectively.  This  model  is  referred  to  as  SBU-4-6  in  the  remainder  of  this  section. 

The  performance  of  model  SBU-4-6  is  compared  to  the  results  of  the  direct  numerical 
simulations  for  the  biaxial,  uniaxial  and  shear  loading  cases.  The  force  displacement  diagrams  in 
addition  to  the  damage  evolution  in  the  interface  and  phase  partitions  are  shown  in  Figs.  3-5.  In 
biaxial  loading,  the  failure  along  the  interface  is  uniform  and  precedes  the  failure  within  the  matrix 
phase.  Upon  interface  debonding,  the  failure  within  the  matrix  propagates  in  the  vertical  and 
lateral  directions.  The  evolution  of  damage  within  the  matrix  partitions  and  the  interface  clearly 
show  that  the  failure  modes  are  accurately  captured  by  SBU-4-6.  A  similar  trend  is  observed  in 
model  predictions  when  subjected  to  uniaxial  (Fig.  5)  and  shear  (Fig.  4)  loading  conditions.  The 
failure  mechanisms  are  captured  with  good  accuracy  when  compared  to  the  reference  direct 
numerical  simulations  of  the  RVE. 


Figure  3.  Stress-strain  Response  and  Damage  Evolution  within  the  RVE  when  Subjected  to 

Uniform  Biaxial  Loading 
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Figure  4.  Stress-strain  Response  and  Damage  Evolution  within  the  RVE  for  Shear  Loading 


Figure  5.  Stress-strain  Response  and  Damage  Evolution  within  the  RVE  when  Subjected  to 
Uniaxial  Tensile  Loading  in  the  Lateral  Direction 


Figure  6  illustrates  the  capability  of  the  proposed  reduced-order  model  in  eliminating  the 
spurious  residual  stresses  in  the  post-failure  regime.  The  spurious  residual  stresses  present  due  to 
the  modeling  errors  associated  with  reduction  of  the  model  order  typically  pollute  the  post-failure 
stress  fields  in  the  macroscopic  analyses,  since  this  effect  partially  constrains  stress  redistribution. 
Figure  6  compares  the  predictions  of  the  proposed  model  along  with  the  predictions  of  an 
eigendeformation-based  homogenization  model  (EFIM  (0+1)  point  model)  proposed  in  Ref.  [45] 
along  with  the  direct  numerical  simulations  when  subjected  to  uniform  biaxial  tension.  The 
matrix-reinforcement  interface  is  assumed  to  remain  continuously  bonded  throughout  the 
simulation.  A  1 -partition  reduced  order  model,  SBU-0-1,  is  adopted.  The  predictions  of  the  EFIM 
(0+1)  point  model  clearly  demonstrate  a  residual  strength  upon  failure  of  the  matrix  partition, 
while  SBU-0-1  eliminates  the  spurious  residual  stresses. 
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Figure  6.  Stress-strain  Response  when  Subjected  to  Uniform  Biaxial  Tensile  Loading 

The  values  of  the  model  parameters  used  in  the  reduced  order  model  are  different  than 
those  of  the  direct  numerical  simulations.  The  objective  of  the  proposed  reduced  order  model  is  to 
capture  the  failure  mechanisms  within  the  heterogeneous  material  in  a  computationally  efficient 
and  accurate  manner.  The  simulations  conducted  in  this  section  demonstrate  that  the  main  failure 
mechanisms  are  captured  with  reasonable  accuracy  with  the  reduced  order  model.  The  model 
parameters  for  the  reduced  order  model  are  computed  by  minimizing  the  discrepancy  between  the 
ultimate  strength  predicted  by  the  proposed  reduced  order  model  and  the  direct  numerical 
simulations  in  a  least  square  sense.  From  the  validation  perspective,  the  reduced  order  model  can 
be  adopted  to  predict  the  response  of  heterogeneous  systems  by  calibrating  the  material  parameters 
of  the  damage  models  directly,  based  on  experimental  observations.  A  general  discussion  and 
methodologies  for  calibration  and  validation  procedures  for  multiscale  models  are  discussed  in 
Ref.  [46], 

Crack  Propagation  in  a  Beam  Subjected  to  Three-Point  Bending 

We  consider  a  three-point  bending  of  a  cracked  composite  plate.  The  predictions  of  the 
SBU-4-6  model  are  compared  to  a  fine  scale  finite  element  model,  which  consists  of  256  RVEs 
described  in  previously  in  this  section.  The  macroscale  mesh  for  the  multiscale  model  consists  of 
256  4-noded  quadrilaterals.  The  volume  fraction  of  the  circular  inclusions  is  30%.  The  circular 
inclusions  are  assumed  to  be  isotropic  and  linear  elastic  with  E  =  200  GPa  and  v  =  0.3  .  Damage 
processes  are  considered  within  the  central  third,  and  the  matrix  material  is  assumed  to  be  linear 
elastic  in  the  remainder  of  the  plate.  The  elastic  properties  of  the  matrix  material  are  E  =  60  GPa 
and  v  =  0.3  .  The  initial  vertical  matrix  crack  is  assumed  to  extend  l/8th  of  the  plate  width. 

Figures  7a  and  7b  illustrate  the  propagation  of  the  initial  matrix  crack  and  damage  within 
the  matrix  as  predicted  by  the  direct  numerical  simulation  and  the  SBU-4-6  model.  In  these 
simulations,  the  interface  between  the  matrix  and  the  inclusions  are  assumed  to  remain  fully 
bonded  for  the  duration  of  the  loading.  The  propagation  of  the  initial  crack  is  arrested 
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approximately  halfway  through  the  plate  thickness  when  shear  cracks  develop  at  the  edges  of  the 
applied  loading.  Figure  7a  shows  the  damage  state  within  the  third  phase  partition  depicted  in  Fig. 
2b.  A  comparison  of  the  reaction  force-applied  displacement  curves  of  the  numerical  simulation 
and  the  proposed  multiscale  model  is  shown  in  Fig.  8.  The  reduced-order  model  slightly 
over-predicts  the  strength  of  the  composite  plate.  The  errors  associated  with  the  SBU-4-6  are  due 
to  the  blunting  of  the  response  fields  across  the  failure  paths  within  the  microstructure  by  the 
model-reduction  methodology.  SBU-4-6  successfully  captures  the  propagation  and  arrest  of  the 
initial  crack  and  subsequent  shear  crack  formation  with  reasonable  accuracy. 
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Figure  7.  Damage  Profile  of  the  3-point  Bending  Beam  Specimen  at  the  Onset  of  Shear 
Fracture  in  the  Absence  of  Interface  Debonding  Effects:  (a)  The  Prediction  of  the  SBU-4-6 
Model;  (b)  Damage  Distribution  Predicted  by  the  Direct  Numerical  Simulation 
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Figure  8.  Comparison  of  the  Load-deflection  Curve  between  the  Multiscale  SBU-4-6  Model 
and  the  Direct  Numerical  Simulation  in  the  Absence  of  Interface  Debonding  Effects 


Figures  10a  and  10b  show  the  failure  of  the  three-point  bending  plate  in  the  presence  of 
interface  effects.  The  direct  numerical  simulation  with  the  fine  mesh  shows  that  the  path  of  crack 
propagation  is  significantly  altered  when  the  inclusion-matrix  debonding  is  considered.  The  path 
of  crack  propagation  displays  a  more  jagged  pattern  with  interaction  between  the  matrix  and 
interface  cracks.  Figure  10a  displays  the  state  of  interface  damage  across  the  macroscale.  The 
extent  of  interface  damage  is  predicted  by  the  SBU-4-6  model  with  reasonable  accuracy.  The 
comparison  of  the  applied  force-deflection  curve  predicted  with  the  proposed  multiscale  model 
and  the  direct  numerical  simulations  is  shown  in  Fig.  9.  The  degradation  effect  of  interface 
debonding  on  the  overall  performance  of  the  plate  is  predicted  by  the  proposed  multiscale  model 
with  good  accuracy. 
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Figure  9:  Comparison  of  the  Load-deflection  Curve  between  the  Multiscale  SBU-4-6  Model 
and  the  Direct  Numerical  Simulation  in  the  Presence  of  Interface  Debonding  Effects 
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Figure  10.  The  Deformed  Configuration  of  the  3-point  Bending  Beam  Specimen  in  the 
Presence  of  Interface  Debonding  Effects:  (a)  The  Prediction  of  the  SBU-4-6  Model;  (b) 
Crack  Profile  Predicted  by  the  Direct  Numerical  Simulation 

It  is  important  to  ascertain  the  reduced  order  model’s  computational  performance,  but  one 
should  note  that  a  running  time  comparison  between  the  reduced  order  model  and  the  direct 
numerical  simulation  depends  highly  on  the  specific  features  of  the  problem.  In  this  case,  the  direct 
numerical  simulation  had  a  mesh  containing  approximately  90,000  elements.  With  this 
configuration,  the  direct  numerical  simulation  ran  over  several  hours  in  comparison  to  the 
SBU-4-6  model  which  ran  in  several  minutes.  If  the  RVEs  had  more  complex  geometry  requiring 
higher  mesh  density,  the  direct  numerical  simulation  would  have  had  many  more  elements  and 
thus  a  slower  running  time.  However,  if  the  same  number  of  failure  paths  were  used  in 
constructing  a  reduced-order  model  for  the  more  complex  microstructural  geometry,  the  running 
time  of  the  reduced  order  model  would  remain  more  or  less  the  same.  In  general,  increasing  the 
number  of  elements  required  in  meshing  a  single  RVE  increases  the  relative  performance  of  the 
reduced  order  model.  Though  even  here,  with  such  a  simple  RVE,  the  performance  was  increased 
by  an  order  of  magnitude. 

3.3  Multiple  Spatio-Temporal  Scale  Modeling  Of  Composites  Subjected  To  Cyclic  Loading 

A  plethora  of  experimental  investigations  in  the  past  few  decades  have  shed  light  into  the 
failure  mechanisms  in  fiber  reinforced  composites  subjected  to  cyclic  loading  (e.g.,  [16,  31], 
among  many  others).  From  the  modeling  perspective,  continuum  damage  mechanics  and  fracture 
mechanics  models  are  typically  employed  to  describe  failure  under  cyclic  loading.  Fracture 
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mechanics  based  approaches  rely  on  incorporation  of  distinct  cracks  at  the  scale  of  the  structure 
[40],  at  the  scale  of  the  constituents  [49],  or  both  [9].  A  fracture  mechanics  based  crack 
propagation  criteria  (e.g.,  Paris  law)  and  a  numerical  methodology  for  crack  propagation  such  as 
mesh  refinement  [34],  virtual  crack  closure  [36],  cohesive  zone  [37],  or  the  extended  finite  element 
method  [35]  are  employed  to  describe  fracture  events  within  the  composite  material.  In  the 
continuum  damage  mechanics  approach,  failure  is  described  as  the  initiation  and  growth  of  diffuse 
damage  (e.g.,  microcrack  density)  typically  represented  using  internal  state  variables.  The 
evolution  of  diffuse  damage  as  a  function  of  loading  history  is  modeled  within  the  nonlinear, 
path-dependent  constitutive  modeling  context  by  employing  micromechanically-informed 
damage  evolution  models  such  as  the  critical  element  model  [50,  59]  and  others  [2,  52,  54]. 

Modeling  complex  failure  mechanisms  and  their  interactions  in  composite  structures 
subjected  to  cyclic  loading  is  a  multiscale  problem  in  space  and  time.  Multiple  spatial  scales  exist 
since  many  failure  mechanisms  initiate  and  grow  at  the  scale  of  the  composite  constituents  defined 
by  the  representative  volume  of  the  composite,  whereas  the  overall  failure  is  assessed  at  the  scale 
of  the  structure  or  structural  component.  Multiple  temporal  scales  exist  because  of  the  disparity 
between  the  characteristic  loading  period,  which  may  be  on  the  order  of  seconds,  and  the  overall 
life  of  the  structure  which  may  be  on  the  order  of  years.  The  computational  homogenization 
method  [32,  61]  based  on  mathematical  homogenization  theory  [6,  10,  51,  57]  is  a  powerful 
multiscale  modeling  approach,  which  has  been  applied  to  nonlinear  solid  mechanics  problems 
involving  multiple  spatial  scales  including  the  failure  of  composite  materials. 

Straightforward  application  of  the  computational  homogenization-based  modeling  to 
evaluate  the  cyclic  response  of  composite  structures  is  prohibitive  due  to  the  tremendous 
computational  cost  associated  with  solving  a  two-scale  nonlinear  problem  in  space  for  large 
number  of  time  steps  necessary  to  evaluate  life  under  cyclic  loading.  This  difficulty  is  addressed 
using  two  approaches:  (a)  by  introducing  reduced  order  (meso-mechanical)  models  that  can 
represent  the  small-scale  response  at  a  fraction  of  the  cost  without  significantly  compromising  the 
solution  accuracy;  and,  (b)  by  introducing  cycle-stepping  methodologies  that  eliminate  the  need  to 
resolve  each  load  cycle  throughout  the  life  of  the  structure  or  structural  component.  Reduced-order 
models  based  on  transformation  field  analysis  [19],  proper  orthogonal  decomposition  [63], 
eigenstrains  [23],  and  others  [1,  29,  39]  have  brought  significant  progress  to  reduced-order 
modeling  in  the  presence  of  multiple  spatial  scales.  The  previous  section  along  with  the  recent 
work  of  Oskay  and  coworkers  [45,  17]  propose  a  reduced  order  computational  homogenization 
framework  based  on  the  eigendeformation  idea  which  provides  (a)  the  ability  to  model  multiple 
failure  mechanisms  at  the  micro  structure  including  matrix  and  fiber  cracking,  and  interfacial 
debonding;  and,  (b)  a  hierarchy  of  reduced  order  models  that  can  be  adapted  to  meet  accuracy 
needs.  The  tyranny  of  temporal  scales  is  addressed  by  employing  cycle-jump  technique  [47]  or 
computational  homogenization-based  temporal  multiscale  modeling  [43,  44,  22]  that  has  been 
employed  in  the  context  of  single  spatial  scale  continuum  damage  mechanics  and 
damage-plasticity  models.  More  recently,  Fish  and  coworkers  applied  cycle-jump  techniques  to 
investigate  the  fatigue  life  of  composites  [25,  26],  Despite  progress,  computational  multiple 
spatio-temporal  scale  modeling  for  accurate,  efficient  and  reliable  prediction  of  failure  in 
composite  structures  subjected  to  cyclic  loading  conditions  remains  to  be  a  challenge. 

In  this  section,  a  new  multiple  spatio-temporal  scale  model  for  prediction  of  cyclic  failure 
in  composite  materials  is  presented.  The  capabilities  of  the  proposed  model  are  demonstrated 
using  a  suite  of  experiments  conducted  on  graphite  fiber  reinforced  epoxy  composites.  The 
proposed  model  is  devised  using  the  computational  homogenization  theory  with  multiple  spatial 
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and  temporal  scales.  The  idea  of  almost  periodicity  of  the  response  fields  at  the  temporal  scales 
[43]  is  employed  to  account  for  the  presence  of  irreversible  damage  fields  that  violate  the 
commonly  assumed  periodicity  conditions  of  the  response  fields.  The  proposed  model  employs  a 
reduced  order  modeling  approach  at  the  spatial  domain  using  the  eigendeformation  based 
homogenization  with  symmetric  coefficients,  and  an  adaptive  time  stepping  strategy  based  on  the 
modified  multistep  method  to  efficiently  evaluate  the  response  of  a  structural  component  by 
resolving  only  a  fraction  of  the  total  number  of  cycles  to  failure.  The  proposed  multiscale  model  is 
calibrated  based  on  a  suite  of  experiments  conducted  on  graphite  fiber  reinforced  epoxy 
(IM7/977-3)  composite  specimens  under  monotonic  and  cyclic  loads,  and  validated  against  an 
independent  set  of  experiments.  The  novel  contributions  found  in  this  section  are  two-fold:  To  the 
best  of  the  authors’  knowledge,  this  is  the  first  attempt  to  concurrently  employ  computational 
homogenization  method  with  multiple  temporal  and  spatial  scales  for  failure  modeling  of 
heterogeneous  materials  subjected  to  cyclic  loading.  The  second-order  adaptive  time  stepping 
methodology  proposed  for  adaptive  error  control  in  time  provides  improvements  in  accuracy 
compared  to  the  first  order  time  stepping  approaches  commonly  employed  in  the  context  of  cycle 
jump  and  temporal  homogenization. 

The  remainder  of  this  section  is  organized  as  follows:  Section  3.3.1  describes  the  multiple 
spatio-temporal  problem  setting.  In  Section  3.3.2,  the  computational  approach  used  modeling 
cyclic  failure  behavior  of  composites  is  formulated.  Section  3.3.3  and  3.3.4  provide  the 
implementation  details  and  the  verification  of  the  proposed  modeling  approach,  respectively. 
Section  3.3.5  describes  the  experiments  for  calibration  and  validation  of  the  computational 
approach,  the  calibration  procedure  of  the  material  parameters  based  on  the  experimental  data,  and 
the  validation  of  the  model. 

3.3.1  Problem  Statement 


time 


loadine  history 


structure 


y  =  */<; 


►v,  representative  volume 


load  cycle 


Figure  11.  Multiple  Spatial  and  Temporal  Scales 


We  consider  the  progressive  failure  of  a  composite  structure  subjected  to  cyclic  loading 
conditions.  Let  QczRl/  be  the  domain  of  a  heterogeneous  body,  where  d  —  1,2  or  3  denotes 
the  number  of  space  dimensions.  Q  is  composed  of  the  repetition  of  a  small  periodic 
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representative  volume  element  (RVE),  0  cz  ,  composed  of  two  or  more  distinct  constituent 
materials  as  illustrated  in  Fig.  11.  The  governing  equations  describing  the  failure  of  the 


heterogeneous  body  are  defined  as  ( x  e  Q  and  t  e  [0  ,tf\ ): 

V-rV'7(x,t)  +  b-  (x)=  0  (69) 

aCn  (x,  t)  =  [l  -  orn  (x,  f )] 1/  (x) :  e*n  (x,  t )  =  if  (x) :  [s^  (x,  t)-  nQn  (x,  t)]  (70) 

sin  (x,  t)  =  V'u^  (x,  t)  (71) 

a®  (x,  t)  =  f(n  (cA'7  ,s(\  s°! )  (72) 


where,  u^';  denotes  displacement  field;  cr-'7  the  Cauchy  stress;  £in  the  total  strain; 
ah'7  e  [0,1)  the  scalar  damage  variable;  /A7  =  or'1  £'n  the  inelastic  strain  tensor;  b-  the  body 
force;  and  if  the  tensor  of  elastic  moduli  obeying  the  conditions  of  symmetry  and  positivity.  The 
evolution  of  or'1  is  typically  nonlinear  and  history-dependent,  and  is  provided  in  the  functional 
form  as  a  function  of  strain,  stress  and  additional  state  variables,  s,;/ .  A  superposed  dot  denotes 
the  material  time  derivative;  x  the  position  vector  parameterizing  the  domain  of  the  structure;  t 
the  time  coordinate;  and,  V-(-),  V(-)  and  Vs(-)  the  divergence,  gradient  and  symmetric 
gradient  operators,  respectively.  The  boundary  conditions  prescribed  on  the  body  consist  of  slowly 
varying  and  oscillating  components  as  illustrated  in  Fig.  11. 

uf'7(x,i)  =  u'7(x,i);  xeT ;  fe[0,fo]  (73) 

•n  =  t'7(x,t);  x e T,;  te[0,to]  (74) 

where,  u  and  t  are  the  prescribed  displacements  and  tractions  on  the  boundaries  ri(  and  T, , 
respectively  ( T  =  u  and  T1(  nT,  =  0);  and,  n  is  the  unit  normal  to  T, .  The  period  of 
oscillations  is  taken  to  be  slow  enough  that  inertial  forces  are  insignificant  and  the  response 
remains  quasi-static.  The  superscripts  C,  and  //  indicate  that  the  response  fields  fluctuate  in 
space  and  time,  respectively.  Double  superscript  indicates  a  response  field  that  fluctuates  in  both 
time  and  space.  The  spatial  fluctuations  arise  due  to  the  fluctuating  material  properties  within  the 
RVE,  whereas  temporal  fluctuations  are  due  to  the  fast  oscillatory  component  of  the  loading.  The 
fluctuating  spatio-temporal  response  is  represented  by  introducing  microscopic  and 
microchronological  scales  parameterized  by  y  =  x/<^  and  r  =  t/// ,  respectively;  and,  o<c  =  i 
and  0  <  r]  =  1  are  scaling  parameters.  The  original  response  fields  that  fluctuate  in  space  and  time 
are  expressed  as: 

<f>in(x,t)=<f>(x,y(x),t,T(t))  (75) 

where,  (j)  denotes  an  arbitrary  response  field.  The  macroscopic  spatial  derivatives  of  a  response 
field  are  obtained  through  the  chain  rule: 

V^(x)  =  VI0(x,y)+-Lv  0(x,y)  (76) 

in  which,  Vx(-)  and  V  (•)  are  gradient  operators  with  respect  to  macroscopic  and  microscopic 

coordinates,  respectively.  All  response  fields  are  assumed  to  be  locally  periodic  with  respect  to  the 
microscopic  coordinates  within  the  RVE  throughout  the  defonnation  process: 

<j){x,  y )  =  (j){x ,  y  +  ky ) ,  where,  y  denotes  the  periods  of  the  microstructure;  and  k  is  a  d  x  d 
diagonal  matrix  with  integer  components.  We  consider  the  following  spatial  homogenization 
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operator: 


*>=(*),=  (77) 
where,  1 0 1  denotes  the  volume  of  the  RVE. 

The  macrochronological  derivative  of  a  response  field  is  expressed  using  the  chain  rule  as: 

f  (0  =  4{t,  t)  =  fa  (t ,t)+  -  $  (t,  t)  (78) 

n 

where,  a  comma  followed  by  a  subscript  variable  t  and  r  denotes  the  partial  time  derivative 
with  respect  to  the  macrochronological  and  microchronological  coordinates,  respectively.  In 
contrast  to  spatial  variability,  local  periodicity  is  not  a  valid  assumption  for  the  response  fields  that 
vary  in  time.  This  is  due  to  the  presence  of  irreversible  mechanisms  associated  with  damage 
accumulation  during  a  load  cycle.  The  response  fields  are  therefore  assumed  to  be  almost  periodic, 
which  implies  that  at  neighboring  points  in  a  temporal  domain  homologous  by  the  load  period,  the 
change  in  the  value  of  a  response  function  is  small  but  does  not  vanish  [44,  22],  Let: 

(</>)T  =— Jor>(x,y,t,r)rfr  (79) 

ro 

denote  the  temporal  averaging  operator;  and,  r0  denote  the  period  of  scaled  cyclic  load  ( 
r  e  [0,  r0] ).  In  the  rate  form,  the  almost  periodic  temporal  homogenization  operator  is  [43]: 

f (0  =  M0),r  =  fit  (0  +  ‘ 1>ap  (0  (8°) 

which  satisfies  the  weak  convergence  property  with  respect  to  an  arbitrary  homogenization 
operator,  ( ■) ;  and,  (f)ap  =  ( j)T)T / 77 .  Following  the  ideas  of  spatial  homogenization  theory,  the 

natural  choice  for  the  temporal  homogenization  operator  is  the  temporal  averaging  operator 
provided  in  Eq.  79.  From  the  computational  perspective,  it  is  more  convenient  to  choose  a 

fixed-point  operator  that  has  the  distributive  property  (i.e.,  (j)  =  t//C  ->(/)=  ij/£)  as  evidenced  by 
the  ensuing  formulation.  In  this  study,  ^(x, y,  t)  =  (f)(\, y, /,())  is  adopted. 

3.3.2  Computational  Model 

This  section  describes  the  multiscale  spatio-temporal  modeling  approach  to  evaluate  the 
failure  response  of  heterogeneous  bodies  subjected  to  cyclic  loading  governed  by  Eqs.  69-74.  The 
eigendeformation-based  homogenization  method  with  symmetric  coefficients  [17]  is  employed  to 
address  the  multiple  spatial  scales,  whereas  temporal  homogenization  with  almost  periodic  fields 
is  employed  to  efficiently  predict  the  life  of  a  structure  without  resorting  to  full  cycle-by-cycle 
analysis.  The  cyclic  damage  evolution  law  employed  to  idealize  the  failure  response  of  composite 
constituents  is  presented. 

Multiple  Scale  Model 

We  start  by  expressing  the  displacement  field  of  the  heterogeneous  body  using  a  two-scale 
asymptotic  expansion: 

u(x,y,f,r)=  u(x,t,r)  +  4'u1(x,y,t,r)  (81) 
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in  which,  U  and  u  are  the  macroscopic  and  microscopic  displacement  fields,  respectively.  The 
governing  equations  are  decomposed  into  macroscopic  and  microscopic  problems  through 
asymptotic  analysis  of  the  governing  equations,  which  consists  of  substituting  Eq.  81  into  Eqs. 
69-74,  and  collecting  the  same  order  terms  of  the  resulting  decompositions.  The  two  leading  order 


equilibrium  equations  are: 

0(C‘):  Vy-a(x,y,t,r)=0  (82) 

0(1) :  Vx-o-(x,y,t,r)  +  Vy-o-1(x,y,t,r)+b(x,y)=0  (83) 

where,  Vx  •(•)  and  Vy  ■(■)  are  divergence  operators  with  respect  to  macroscopic  and 

microscopic  spatial  coordinates,  respectively;  and,  the  body  forces  are  assumed  to  remain  constant 
in  time  for  simplicity.  The  first  and  second  order  stress  fields  are: 

cr(x,  y ,  t,  r)  =  [l  -  co(x,  y,  t ,  r)]L(y ) :  [e (x,  t,  r) + Vy  (u1  )(x,  y,  t,  r)]  (84) 

a1  (x,  y ,  t,  r)  =  [l  -  oj(x,  y ,  t,  r)]L(y ) :  [vx  (u 1  )(x,  y ,  t ,  r)]  (85) 


in  which,  Vx(-)  and  V s ( ■ )  are  symmetric  gradient  operators  with  respect  to  the  macroscopic  and 

microscopic  spatial  coordinates,  respectively;  and,  s  =  Vx(u)  is  the  macroscopic  strain  tensor. 
Applying  Eq.  78  to  the  damage  evolution  equations  and  collecting  the  same  order  terms 

yield: 

0(jil):o)T=f0{(7,£,  s)  (86) 

O(\):cot  =f1(a,£,s)  (87) 

where,  the  evolution  functions  f°  and  /'  are  derived  based  on  the  prescribed  evolution  law  f 
The  boundary  data  applied  to  the  heterogeneous  body  is  composed  of  a  slowly  varying  and 
a  periodic  oscillatory  component: 

u,;(x,t)  =  u°(x,t)  +  u1(x,  r)  (88) 

t'?(x,i)  =  t°(x,i)+t1(x,r)  (89) 

Applying  spatial  averaging  (Eq.  77)  to  the  0(1)  equilibrium  equation  (Eq.  83),  exploiting 
the  local  periodicity  of  the  stress  fields  in  space,  and  applying  the  almost  periodic  temporal 
homogenization  operator  yields  the  macroscale  equilibrium  equation.  The  resulting  equilibrium 
equation  along  with  the  spatio-temporally  homogenized  first  order  stress  field,  and  the  boundary 
conditions  provide  the  macrochronological-macroscopic  boundary  value  problem. 

Macrochronological  -  Macroscopic  Problem:  Given:  average  body  force,  b, 
boundary  data  u°  and  t° ,  and  the  solution  of  the  macrochronological-microscopic  problem;  Find 
the  macroscopic  displacement  field,  u,  such  that  (t  s  [0  ,tf] ): 


Vx -cr(x,t)+b(x)  =  0;  xeD 

(90) 

^(x,  t )  =  ((l  -  <5)L(y) :  \s  +  Vy  (u1)^ ;  x  e  Q 

(91) 

n  =  u°(x,0;  xer„ 

(92) 

a(x,t)  -n  =  t°(x,f);  xeT, 

(93) 

Applying  the  temporal  homogenization  operator  to  the  0(£  1 )  equilibrium  equations  (Eq.  82) 
along  with  the  constitutive  equation  for  the  leading  order  stress  (Eq.  84),  and  employing  periodic 


32 

Approved  for  public  release;  distribution  unlimited 


boundary  conditions  for  the  microscale  displacement  field,  the  macrochronological-microscopic 
boundary  value  problem  is  obtained. 

Macrochronological-Microscopic  Problem:  At  a  fixed  macroscale  material  point 
x  e  Q  ;  Given:  macroscale  strain,  s ,  and  the  tensor  of  elastic  moduli,  L  ;  Find  the  displacement 
field,  u1 ,  such  that: 

V,-|l-S)L(y):[?  +  V;(u,)}=0;  y  *9  (94) 

a(x,t)  =  f{v,e,s)+(oap{%ty,  ye0  (95) 

u1  periodicony  ed0  (96) 

The  macro-  and  microscopic  problems  associated  with  the  fast  time  scale  at  a  fixed  slow  time 
coordinate,  t ,  are  obtained  based  on  similar  algebra,  but  without  applying  the  temporal 
homogenization  operator  to  the  governing  equations  and  considering  the  damage  evolution 
equation  with  respect  to  the  fast  time  scale  (i.e.,  Eq.  86).  The  resulting  microchronological  - 
macroscopic  and  microchronological  -  microscopic  problems  are  stated  as  follows: 

Microchronological  -  Macroscopic  Problem:  At  a  fixed  macrochronological  time 

t  g  [0 ,tf\.  Given:  average  body  force,  b,  boundary  data,  u  and  t ,  and  the  solution  of  the 
microchronological-microscopic  problem;  Find  the  macroscopic  displacement  field,  u ,  such  that 
(r  G[0,r0]): 


Vx -cr(x,r,r)+b(x)=  0;  xgQ 

(97) 

:((l-ru)L(y):[^  +  V;(u')])Y;  xgQ 

(98) 

u  =  u(x,f,r);  xeT 

(99) 

crn  =  t(x,F,r);  xgT, 

(100) 

Microchronological-Microscopic  Problem:  At  a  fixed  macroscale  material  point 
x  e  Q  and  a  fixed  macrochronological  time  t  g  [0 ,G];  Given:  macroscopic  strain,  e ,  and  the 

tensor  of  elastic  moduli,  L  ;  Find  the  microscopic  displacement  field,  u1 ,  such  that: 

Vy  '{(l-(y)L(y):  [f  +  Vy(u')]}=  0;  yet?  (101) 

cur(x,y,Gr)  =  /°(cr,^,s);  y  e  0  (102) 

u1  periodicony  g  dO  (103) 

The  macrochronological  and  microchronological  problems  are  coupled  through  the  almost 
periodic  rate  operator  that  defines  the  evolution  of  the  temporally  homogenized  response  fields 
(i.e.,  Eq.  95).  Therefore,  the  evolution  of  the  macrochronological  fields  at  each 
macrochronological  time  coordinate  requires  the  solution  of  the  microchronological  problem 
associated  with  that  time  coordinate.  The  macroscale  problems  at  the  macrochronological  and 
microchronological  time  scales  are  coupled  with  the  respective  microscale  problems  through  the 
constitutive  relationship  (Eq.  98).  The  evaluation  of  the  macroscopic  stress  at  each  macroscopic 
material  point  requires  the  solution  of  the  microscopic  RVE  problem  associated  with  that  material 
point.  When  the  finite  element  method  is  employed  to  evaluate  the  macroscale  problem,  a 
nonlinear  microscale  problem  must  be  evaluated  to  update  the  stress  at  each  integration  point  for 
each  increment  and  iteration  of  every  time  step  of  the  loading  history.  This  is  a  significant 
computational  burden. 
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Reduced  Order  Spatial  Homogenization 


We  employ  the  eigendeformation-based  reduced  order  homogenization  method  with 
symmetric  coefficients  (sEHM)  to  reduce  the  computational  cost  associated  with  evaluating  the 
coupled  nonlinear  micro-  and  macroscopic  problems.  sEHM  is  introduced  in  Section  2  and  a  brief 
summary  is  provided  herein.  The  premise  of  sEHM  is  to  devise  a  low-cost  approximation  to  the 
nonlinear  microscale  boundary  value  problem  defined  over  the  representative  volume  elements 
based  on  the  idea  of  precomputing  certain  microstructural  information  (e.g.,  concentration  tensors, 
localization  operators,  influence  functions)  through  linear  elastic  simulations  prior  to  the  analysis 
of  the  macroscale  structure. 

The  microscale  displacement  field  is  expressed  as: 

u1  (x,  y ,  t,  t)  =  H(y ) :  e  (x,  t,  t  )  ■ +  £  h(y ,  y ) :  //(x,  y ,  t,  r)dy  (104) 

in  which,  H  is  the  elastic  influence  function  (a  third-order  tensor)  obtained  by  substituting  Eq. 

104  into  the  microscale  problem,  and  evaluating  the  microscale  problem  in  the  absence  of  damage; 
and  h  is  the  phase  damage  induced  influence  function  provided  by  the  particular  solutions  to  the 
RVE  problems  obtained  by  substituting  Eq.  104  into  the  microscale  problem,  and  solving  the 
microscale  problem  in  the  presence  of  phase  damage  (i.e.,  // ).  The  governing  equations  and  the 
discrete  approximations  of  the  elastic  and  phase  damage  induced  influence  functions  are  provided 
in  Ref.  [45],  Meso-mechanical  shape  functions  are  employed  to  discretize  damage  and  damage 
induced  inelastic  strain  fields: 


a>,  //}(x,  y,  t)  =  2>w(y){«(r),  //r)}(x,  t,  r ) 


(105) 


r= i 


in  which,  N{yy)  are  the  phase  shape  functions.  We  consider  a  partitioning  of  the  RVE  domain  into 
n  non-overlapping  subdomains,  6yy) ,  such  that  G(y)  n6{A)  =0  if  y  ^  A .  The  phase  shape 
functions  are  taken  to  be  piecewise  constant  functions  forming  a  partition  of  unity  within  the  RVE: 

\\ify  e0(r) 

J  3  (106) 
[0  elsewhere 

and,  oJ'y)  and  ju^  are  damage  variable  and  inelastic  strains  averaged  over  the  partition,  0(y) . 
Substituting  Eqs.  104-106  into  the  microscale  problem  (Eqs.  3.1),  and  considering  variational 
arguments  (see  Ref.  [17]),  the  governing  equation  of  the  microscale  problem  is  reduced  to  the 
following  algebraic  form: 


A(orA):£  + 


YB(aAy) :  //<7) 


7= 1 


=  0  Va  =  l,2,...,H 


(107) 


in  which,  A(aA)  and  B<QAr)  are  coefficient  tensors  computed  as  a  function  of  the  influence 
functions,  H  and  h  as  well  as  the  elastic  properties,  L  .  The  expressions  for  A("A)  and  B(aAr) 
are  found  in  Eqs.  35  and  37  where  (A =  and  (B(aAr))..w  =  F^Ay)  in  this  case.  The 


macroscopic  stress  tensor  is  expressed  in  terms  of  the  partition  average  damage  variable  and 
inelastic  strains  as: 

a  =  2J  -  co(A)  ][l(A)  :  £+  ^p<  :  ^ 


A=1 


a= 1 


(108) 
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in  which,  L(A)  and  P<Aff)  are  coefficient  tensors  whose  expressions  are  found  in  Eqs.  49-50.  The 
macrochronological  counterparts  of  the  reduced  order  microscopic  equilibrium  and  macroscopic 
stress  are  obtained  by  applying  the  almost  periodic  temporal  homogenization  operator  to  Eqs.  107 
and  108. 

Cyclic  Damage  Model 


Strain  Strain 

Figure  12.  Effect  of  the  Cyclic  Damage  Sensitivity  Parameter,  p  ,  on  the  Cyclic  Stress-strain 

Relationship 


Continuum  damage  mechanics  (CDM)  is  employed  to  describe  the  evolution  of  damage 
within  a  phase  partition  (i.e.,  (o(y) ).  In  contrast  to  monotonic  CDM  models,  the  evolution  law  is 
allowed  to  accumulate  damage  at  subcritical  loading  levels  to  permit  sensitivity  to  cyclic  loading. 
Such  a  model  previously  employed  in  Refs.  [25,  43],  is  adopted  in  this  study: 

a>(/)  =gp  d°.^r)  )<*>(y)  >+  where  0 <g  =  ^<1;  y  =  l,2,...,n  (109) 

du n  co  n 

where,  p  is  the  cyclic  damage  sensitivity  parameter;  denotes  MacCauley  brackets;  ®  the 

damage  evolution  law  under  monotonically  increasing  loads;  and  uly)  the  damage  equivalent 
strain: 


v 


(r)  — 


(110) 


in  which,  L(7)  is  the  tensor  of  elastic  moduli  of  the  constituent  occupying  0l  ] ;  and,  siy)  is  the 
average  strain  within  partition,  y .  The  damage  evolution  under  monotonic  conditions  is  idealized 
based  on  a  smooth  evolution  law: 

0^,0) arctan(o(r)u(7) -6(7))+arctan(&<7)) 

y  +  arctan(/f7)) 


where,  a(7)  and  b(y)  are  material  parameters  associated  with  the  constituent  occupying  domain 

0(r) . 

The  cyclic  damage  sensitivity  parameter,  p  ,  which  controls  the  rate  of  damage 
accumulation  with  respect  to  a  load  cycle  as  illustrated  in  Fig.  12,  is  taken  to  be  of  the  form: 

35 

Approved  for  public  release;  distribution  unlimited 


(112) 


(r)  —  J,r)  i  Jr).p)  ,  Sri (,  Sr)  A 

P  (0  ^C1  Gnax  ^  C2  '  Gnax ) 

Cgr) ,  Cg  ] ,  and  c(2r>  are  parameters  that  account  for  the  sensitivity  of  fatigue  strength  to  the 
maximum  applied  loading.  A  quadratic  form  is  chosen  to  model  the  fatigue  behavior  observed  in 
the  experiments.  is  the  maximum  value  of  the  damage  equivalent  strain  during  the  history  of 
the  material  point  as  defined  below. 

^(0  =  max{u(r)(r)|0<r<4  (113) 

T 

Applying  the  chain  rule  in  the  time  domain  (Eq.  78)  to  the  damage  evolution  equation  (Eq. 
109),  collecting  the  terms  based  on  the  order  of  the  temporal  scaling  parameter,  7 ,  and  applying 
the  temporal  homogenization  operator  to  the  0(1)  equation  yields: 

<=/"  =  <114> 

s>T=f'  =  sr 

Microchronological  Problem:  Given:  Coefficient  tensors,  L(A),  P(AQf),  A(qA)  and 
B(aAr) ;  average  body  force,  b,  boundary  data  u  and  t  Find:  at  a  fixed  macrochronological  time 
t  e [0,0]  the  macroscopic  displacement  field,  u,  which  satisfies  (re[0,rj): 

Equilibrium:  Vx  -cf(x,F,z-)+b(x)  =  0;  xeO 
Constitutive  Equation :  a  =  ^[l  -  a>(A)  ][l(a)  :  s+  ^p(Aa) ;  ju(a) 

A=1  a= 1 

jj[l-<y(A)’  AiaA):e  +  YB^aAr):^r)  1  =  0  Va  =  l,2,...,/i 

A=1  [  ]_  7=1  JJ 

Damage  Evolution :  e/;y)  =  f°  =  gp  ^^7)  ~  <  (;.7>  >+ 

Boundary  Conditions  :  u  =  u(x,  t ,  r);  x  e  r,( 
an  =  t(x,f,  r);  xeT, 

3.3.3  Computational  Implementation 

The  microchronological-  and  macrochronological-reduced  order  multiscopic  problems  are 
described  in  the  previous  section.  Given  the  coupling  terms,  the  micro-  and  macrochronological 
problems  are  evaluated  using  the  nonlinear  finite  element  method.  A  commercial  finite  element 
software  (Abaqus)  along  with  the  user  material  subroutine  utility  (UMAT)  is  employed  to  solve 
these  problems.  In  this  section,  we  focus  on  the  implementation  details  of  the  coupling  between  the 
micro-  and  macrochronological  problems.  The  proposed  solution  strategy  is  implemented  with  an 
adaptive  macrochronological  time  stepping  methodology. 

The  macrochronological  system  of  equations  for  evaluation  of  a  macroscale  time  step. 
Macrochronological  Problem:  Given:  coefficient  tensors,  L(A),  P(Aof),  A("A)  and  B(aAr); 
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r/o(gw) 


<uir)  > 


(115) 


average  body  force,  b,  boundary  data,  u°  and  t° ;  the  almost  periodic  damage  function,  coap  . 
Find :  the  displacement  field,  u,  which  satisfies  ( /  e  [0,/  /  ] ): 

Equilibrium:  Vx  •or(x,t)  +  b(x)  =  0;  xeQ 

Constitutive  Equation :  <7  =  ^[l  -  &>(A)  ][l<a)  :  s  +  ^p(Ao,) ;  fi{a) 

A=1  a= 1 


■CO 


(A) 


(aA)  . 


+  Ypia&r)  ://<r) 


Y= 1 


Damage  Evolution:  S(x,t)  =  f  +  coa  =  gp 

duyn 


=  0  Va  =  1,2, 


<D5r)>++®  (x,0; 


a/? ' 


Boundary  conditions:  u  =  u°(x,i);  xeT, 


cr(x,t)-n  =  t°(x,t);  xeT, 

The  overall  solution  strategy  for  the  evaluation  of  the  coupled  multiscale  system  is 
illustrated  in  Fig.  13.  Consider  a  discretization  of  the  macrochronological  time  domain, 

\t0  =0,...,ti_l,ti,ti+1,...,tk  =tf  \  in  which  tf  denotes  the  zth  macrochronological  time  step: 

(116) 

7=1 

where,  At,  =  t,  -tM  and  t0  =  0 .  A  driver  program  (implemented  in  the  Python  programming 

language  for  compatibility  with  Abaqus)  controls  the  execution  of  the  solution  procedure.  At  each 
macrochronological  time  step,  the  time  step  size  as  well  as  the  almost  periodic  damage  field  are 
estimated  and  passed  to  the  macrochronological  -  multiscopic  problem,  which  is  evaluated  for 
macrochronological  response  fields.  The  macrochronological  response  fields  provide  the  initial 
state  of  the  microchronological  -  multiscopic  problem  at  fixed  time  tt  due  to  the  particular  choice 
of  the  temporal  homogenization  operator.  The  point-wise  value  of  the  almost  periodic  damage 
field,  (0^ ,  at  the  current  time  step  is  computed  and  passed  to  the  driver  routine  for  computation  of 
the  next  time  step  size. 
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Figure  13.  Implementation  Strategy  of  the  Coupled  Micro-  and  Macro-chronological 

Problems 


Adaptive  Macrochronological  Time  Stepping 


The  proposed  solution  algorithm  outlined  in  Fig.  13  requires  the  resolution  of  k  cycles 
throughout  the  loading  history,  which  are  evaluated  by  the  microchronological  problem  (Box  13). 
The  accuracy  and  efficiency  of  the  proposed  approach  is  based  on  the  appropriate  selection  of  the 
macrochronological  time  steps,  as  well  as  the  accurate  approximation  of  the  evolution  of  the 
almost  periodic  component  of  the  damage  field.  In  this  study,  the  almost  periodic  damage  fields 
are  approximated  based  on  a  modified  quadratic  multistep  method  [13],  whereas  the 
macrochronological  time  step  size  is  chosen  adaptively  based  on  a  maximum  damage 
accumulation  criterion.  It  is  assumed  that  damage  accumulation  is  primarily  due  to  cyclic  loading 
and  slow  loading  component  remain  smooth  during  the  loading  history. 

Let  D'(t)  be  a  matrix  of  almost  periodic  damage  rate  fields: 


D'(0 


a>. 


(i) 

ap 


>0 


(117) 


in  which,  ng  denotes  the  total  number  of  integration  points  within  the  macroscopic  domain;  x„ 
denotes  the  value  of  the  function  at  integration  point  g .  Considering  a  smooth  damage  growth 
between  macrochronological  steps  i  - 1  and  i  + 1 ,  the  evolution  of  D'(t)  is  approximated  by  a 
linear  function  around  tt : 


D’(t,  +  At)  «  P,-(At)  =  b,  +  c,At;  At  e[-At,,Ati+l]  (118) 

The  coefficients  of  the  linear  approximations  are  obtained  by  imposing  the  following  conditions 
on  P;. : 

P,(0)  =  D’(t,) 

P,(-At,)  =  D'(tM) 
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(119) 

(120) 


Substituting  Eqs.  1 19,  and  120  into  Eq.  118,  the  rate  of  the  almost  periodic  damage  fields  is 
obtained  as: 

D'(/„,)»P(A/w)  =  D'((,)  +  ^-(D'(/,)-D'((,-l))  (121) 

At,. 

Let  AD.+1  define  the  matrix  of  cyclic  loading  induced  damage  accumulated  between  time 
t,  and  tM  .  Employing  the  linear  approximation  provided  by  Eq.  118  yields: 

ADw«E>,(A()<(A/  =  A»wD'(»,)  +  3(D'ft)-D'ft-1))  (122) 

J0  2Atj 

The  macrochronological  step  size  ( At,+1 )  is  chosen  such  that  the  maximum  damage  accumulation 
within  a  step  does  not  exceed  a  threshold  value: 

||AD,„j„„<A£U  (123) 

where,  A Z)  is  the  damage  accumulation  threshold,  and  ll-ll  denotes  the  matrix  max  norm. 

The  loss  of  load  carrying  capacity  during  a  microchronological  problem  following  a  long 
macrochronological  step  indicates  a  possible  underestimation  of  damage  accumulation  during  the 
macrochronological  time  step.  This  is  due  to  the  deviation  of  the  damage  accumulation  from  the 
piecewise  quadratic  approximation  within  the  time  step.  When  failure  is  observed,  the  previous 
macrochronological  step  size  is  shortened,  and  the  previous  macrochronological  time  step  is 
repeated.  The  time  step  shortening  is  repeated  until  desired  accuracy  of  the  time-to-failure  is 
achieved.  The  overall  algorithm  is  presented  below: 

1 .  Initialize  algorithm:  i  =  0 . 

2.  Evaluate  macrochronological  step  (Box  4). 

3.  Evaluate  microchronological  problem  at  t  =  0  (Box  13). 

4.  Set  i  =  l;  At,  =  r/r0 ,  where  77  r0  is  the  duration  of  a  single  load  cycle;  D'(t,)  =  D'(0) . 

5.  Evaluate  macrochronological  step. 

6.  Evaluate  microchronological  problem  at  t  =  tx. 

7.  Set  i  =  2 

8.  While  7  <  k  : 

(a)  Calculate  At,,  using  Eq.  123. 

(b)  Calculate  D'(t,.)  using  Eq.  121. 

(c)  Evaluate  macrochronological  step. 

(d)  Evaluate  microchronological  problem  at  t  =  t, . 

(e)  If  {failure  event  &  At,./t,  >  to/  }: 

i.  At,.  <—  max(77r0,cAt,.)  and  go  to  step  (b). 

(f)  Else  if  {failure  event  &  A t./t.  <toI  }: 

i.  Structural  failure:  Stop  algorithm. 

(g)  /  <—  7  + 1 

9.  End. 

The  proposed  algorithm  is  initiated  by  evaluating  two  macrochronological  time  steps  and  two 
microchronological  problems  that  corresponds  to  the  first  two  load  cycles,  noting  that  a  single  load 
period  is  denoted  by  77  r0 .  The  almost  periodic  field  is  taken  to  vary  linearly  between  the  first  and 
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second  load  cycles,  and  quadratically  between  the  remaining  macrochronological  time  steps. 

0  <  c  <  1  is  the  step  cutback  factor  when  a  failure  event  is  detected.  The  failure  event  is  defined  as 
a  loss  of  load  carrying  capacity  of  the  structure,  detected  as  lack  of  convergence  that  occurs  during 
the  evaluation  of  the  microchronological  problem.  The  macrochronological  time  step  size  cannot 
be  smaller  than  a  single  load  period  (  =  rjr0 ).  Structural  failure  is  taken  to  occur  when  the  failure 

event  is  detected,  and  the  ratio  of  the  current  macrochronological  time  step  size  and  the  current 
macrochronological  time  is  less  than  a  specified  tolerance  ( tol )  value.  The  cutback  iterations 
when  a  failure  event  is  detected  are  controlled  by  the  Python  driver  program,  which  employs  the 
restart  capability  of  Abaqus  to  iterate  the  evaluation  of  macro-  and  micro-chronological  problems 
until  convergence. 

Improved  Adaptive  Stepping  Criterion 


An  improved  adaptive  step  size  criterion  can  be  developed  by  considering  adaptive 
Runge-Kutta  methods  used  in  solving  initial  value  problems.  Adaptive  Runge-Kutta  methods 
estimate  the  truncation  error  by  comparing  an  order  n  method  to  an  order  n- 1  method  while 
adapting  the  step  size  according  to  this  estimate.  In  similar  manner,  we  consider  a  smooth  damage 
growth  between  macrochronological  steps  i  and  i  + 1  of  one  order  less  than  Eq.  1 18  (i.e.  a 
constant  function). 

D'(f,  +  A*)  *  Pf(A0  =  bf;  At  e  [0,  Ati+1)  (124) 

b,  is  attained  by  imposing  Eq.  1 19  on  Eq.  124. 


We  reach  an  expression  for  P  . 


P,(0)  =  D'(O  =  b, 


(125) 


P,(At)  =  D'(t,)  (126) 

In  parallel  with  Eq.  122,  the  matrix  of  cyclic  loading  induced  damage  accumulated  between  times 
t,  and  tM  is  approximated  by  integrating  Eq.  126  from  0  to  Ati+I . 


/•At. , , 

AD,.+1  «  £  1P/(At)dAt  =  D'(t,  )AtI+1 


(127) 


Let  ADr_j  be  the  second  order  approximation  to  AD.+1  in  Eq.  122  and  AD)+I  be  the  order  one 


approximation  of  AD(+1  in  Eq.  127.  The  step  size  is  adaptively  chosen  so  that  the  difference 
between  ADf+1  and  AD:+1  is  less  or  equal  to  a  tolerance  denoted  tol . 


tol> 


|AD:+1  -  AD: 


+i 


At: 


—  w+i 


2  At, 


|D'«,)-D'«,„)|L 


Assuming  equality  in  Eq.  128,  a  closed  form  expression  is  reached  for  the  step  size  A tx 


+1  • 


Af+,  = 


2  tol  At,. 


(128) 


(129) 


|D'(6)-D'(tM)|max 

To  utilize  the  new  adaptive  stepping  criterion  within  the  previously  stated  algorithm,  replace  Eq. 
123  in  step  8a  with  Eq.  129.  Eq.  129  is  ill-defined  if  |D'(t,-)  -  D'(6-i)|max  =  0 .  Therefore,  if 

|D'(6)  -  D'(iw)|  <  s  ,  Eq.  123  is  solved  instead  to  determine  A tM  .  s  is  a  small  numerical 

parameter. 
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3.3.4  Model  Verification 


The  proposed  adaptive  macrochronological  time  stepping  strategy  is  verified  by  comparing 
the  performance  with  direct  cycle-by-cycle  simulations.  In  the  direct  cycle-by-cycle  analysis,  each 
load  cycle  throughout  the  loading  is  resolved,  without  resorting  to  the  multiple  temporal  scale 
strategy.  The  simulations  were  conducted  on  a  unidirectionally  fiber-reinforced  matrix  unit  cell. 
The  fiber  is  taken  to  remain  elastic  throughout  the  loading  period,  whereas  the  damage 
accumulates  within  the  matrix  as  a  function  of  loading  cycles.  The  failure  response  in  the  matrix 
phase  is  approximated  using  a  4-partition  reduced-order  model.  A  uniform  tensile  strain  is  applied 
transverse  to  the  reinforcement  direction.  The  strain  amplitude  is  varied  between  zero  and  the 
maximum  strain  throughout  the  loading  history. 

Figure  14  illustrates  the  variation  of  the  damage  variables  within  the  unit  cell  as  a  function 
of  applied  loading  cycles.  Only  two  of  the  four  matrix  partitions  show  appreciable  damage.  The 
proposed  model  simulations  are  conducted  by  setting  A Dmax  =  0.5%,  1%,  and  2%.  A  cutback 
factor  of  0.5  is  employed  (c  =  0.5  ).  The  proposed  adaptive  time  stepping  strategy  required  109, 
68,  and  43  resolved  microchronological  load  cycles  for  ADmax  =0.5%,  1%,  and  2%,  respectively, 
compared  to  900  cycles  resolved  in  the  direct  cycle-by-cycle  approach.  The  model  captures  the 
failure  response  with  good  accuracy,  particularly  when  A Dmax  =  0.5%. 


Figure  14.  Comparison  of  the  Cyclic  Damage  Accumulation  Computed  using  the  Direct 
Cycle-by-Cycle  Approach  and  Proposed  Multiscale  Model  with  Adaptive  Time  Stepping 
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Figure  15.  Comparison  of  the  Cyclic  Damage  Computed  using  the  Direct  Cycle-by-Cycle 
Approach  and  Proposed  Model  with  Adaptive  Time  Stepping  with  Smaller  Loading 

Amplitude 


The  efficiency  of  the  proposed  approach  is  further  illustrated  by  conducting  simulations 
when  the  unit  cell  is  subjected  to  a  slightly  smaller  loading  amplitude.  In  this  analysis,  the 
cycle-to-failure  of  the  first  partition  is  677  (in  contrast  to  377  of  the  previous  simulations).  The 
simulations  are  conducted  for  1800  cycles.  Figure  15  illustrates  the  variation  of  damage  variables 
within  the  unit  cell  as  a  function  of  applied  loading  cycles  as  computed  by  the  reference 
simulations  and  the  multiple  spatio-temporal  model  with  adaptive  time  stepping  methodology. 
The  total  resolved  cycles  of  the  proposed  adaptive  time  stepping  strategy  remain  largely  the  same  ( 
A Dmax  =  0.5%,  1%,  and  2%  are  1 18,  72  and  47  respectively),  pointing  to  significant 
computational  advantage  in  high  cycle  failure  conditions. 

Improved  Adaptive  Stepping  Criterion  Verification 

In  Fig.  16,  the  improved  adaptive  stepping  criterion  is  verified  by  comparing  three  different 
values  of  tol  to  a  direct  cycle-by-cycle  simulation  for  the  same  loading  case  used  in  Fig.  14. 
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Figure  16.  Comparison  of  the  Cyclic  Damage  Accumulation  Computed  using  the  Direct 
Cycle-by-Cycle  Approach  and  the  Multiscale  Model  with  Improved  Adaptive  Time 

Stepping 


The  three  values  of  tol  were  chosen  to  produce  similar  accuracy  to  the  three  values  of  A Dmax  in 

Fig.  14  (i.e.  0.5  %,  1%,  and  2  %).  The  selected  tol  values  were  0.001 ,  0.002  ,  and  0.004  , 
respectively.  When  comparing  the  new  adaptive  time  stepping  criterion  with  tol  =  0.001  to 
A Dmax  =  0.5% ,  the  same  accuracy  was  maintained  while  resolving  only  53  loading  cycles  as 

compared  to  resolving  109  loading  cycles  with  the  previous  adaptive  stepping  criterion,  a 
significant  computational  savings. 


3.3.5  Assessment  of  Model  Capabilities 

The  capabilities  of  the  proposed  multiple  spatio-temporal  methodology  are  assessed 
through  the  investigation  of  graphite  fiber-reinforced  epoxy  composites  (i.e.,  IM7/977-3) 
subjected  to  cyclic  loading.  This  section  presents  the  experiments  conducted  to  study  the  cyclic 
response  of  the  composite;  calibration  of  the  model  parameters  based  on  monotonic  and  cyclic 
experiments,  and;  validation  of  model  predictions  based  on  acoustic  emission  testing. 

Experiments 

A  suite  of  experiments  was  conducted  to  calibrate  the  material  parameters  and  assess  the 
validity  of  the  proposed  multiscale  model.  Composite  specimens  with  three  separate  layups  of 
unidirectional  laminae  were  tested  under  uniaxial  monotonic  and  cyclic  loading  conditions:  (a) 
Zero  degree  specimens  consist  of  eight  unidirectional  plies  with  fibers  oriented  parallel  to  the 
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loading  direction;  (b)  Ninety  degree  specimens  consist  of  sixteen  unidirectional  plies  with  the 
fibers  oriented  perpendicular  to  the  loading  direction;  and  (c)  Quasi-isotropic  specimens  with  the 
layup  of  [+45,0,-45,90]2i .  Specimen  configurations  are  summarized  in  Table  1.  The  mean  fiber 

volume  fraction  of  the  specimens  is  65.6%,  which  was  determined  based  on  acid  digestion  testing. 
The  results  of  experiments  conducted  on  zero  and  ninety  degree  specimens  are  employed  in  the 
calibration  of  the  parameters  of  the  proposed  multiscale  model,  whereas  the  quasi-isotropic 
specimens  are  employed  in  the  validation  analyses. 

Table  1.  IM7/977-3  Specimen  Dimensions 


Fiber 

Orientation 

Number  of 
Plies 

Length  [mm] 

Width  [mm] 

Thickness  [mm] 

0° 

8 

250 

13 

1 

90° 

16 

177 

25 

2 

Quasi-iso. 

16 

250 

25 

2 

Figure  17.  Failure  Profiles  when  Subjected  to  Monotonic  Loading,  (a)  Zero-degree 
Specimens;  (b)  Ninety  Degree  Specimens;  (c)  Quasi-isotropic  Specimens 


Acoustic  emission  (AE)  was  used  to  detect  failure  events  within  the  quasi-isotropic  layups. 
In-situ  AE  activity  was  recorded  on  a  Micro-II  Digital  AE  System  produced  by  Physical  Acoustics 
Corporation.  In  the  AE  technique,  the  stress  waves  produced  by  the  sudden  release  of  strain  energy 
during  localized  failure  events  are  identified  and  recorded  as  hits.  Appropriate  signal  conditioning 
parameters  are  identified  based  on  an  AE  calibration  study  prior  to  testing.  A  threshold  wave 
amplitude  of  48  dB  enables  the  separation  of  all  valid  failure  events  from  ambient  noise. 

The  first  set  of  experiments  conducted  is  on  unidirectional  unnotched  tension  specimens 
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tested  under  monotonic  displacement  control  on  an  MTS  universal  testing  machine  according  to 
ASTM  Standard  D3039  [3],  A  grip  pressure  of  500  psi  was  applied  to  prevent  slipping  without 
crushing  the  composite.  One  axial  and  one  transverse  strain  gage  were  mounted  on  each  specimen 
to  determine  Poisson’s  ratio.  A  one-inch  extensometer  was  used  to  accurately  measure  the  axial 
stiffness.  All  monotonic  tests  were  conducted  at  a  constant  displacement  rate  of  1 .27  mm/min. 
Thirteen  zero  degree  specimens,  seventeen  ninety  degree  specimens,  and  seven  quasi-isotropic 
specimens  were  subjected  to  uniaxial  tension  up  to  failure  to  ensure  repeatibility.  A  moderate 
degree  of  modulus  and  strength  scatter  is  observed  in  the  experiments.  Figure  17  illustrates  the 
failure  patterns,  which  show  destructive  fiber  failure  in  the  zero  degree  specimens, 
matrix-dominated  failure  in  the  ninety  degree  specimens,  and  combined  matrix  and  fiber  failure  in 
the  quasi-isotropic  specimens.  All  non-zero  plies  in  the  quasi-isotropic  specimens  showed 
matrix-dominated  failure,  while  the  zero  degree  plies  showed  fiber  failure.  The  elastic  and  strength 
properties  observed  in  the  experiments  are  summarized  in  Table  2.  A  Poisson’s  ratio  of  0.3 16  with 
a  standard  deviation  of  0.039  was  observed  by  placing  a  strain  gage  perpendicular  to  the  loading 
on  the  0°  specimens.  The  shear  modulus,  Gu ,  was  determined  to  be  4.66  MPa  with  a  standard 
deviation  of  0.61  using  additional  tension  experiments  conducted  on  composite  specimens  with  a 
+/-45  °  layup  according  to  the  procedure  described  in  ASTM  D35 18  [5], 


Table  2:  Calibrated  Elastic  Parameters  of  the  Composite  Constituents;  Observed  and 
Simulated  Elastic  Parameters  of  the  Overall  Composite  (Standard  Deviation  in  Parenthesis) 


Layup 

Young’s  Modulus, 

E  [GPa] 

Failure  Strength, 
crf  [MPa] 

0° 

158(13) 

2,841  (296) 

90° 

8.644  (0.712) 

63  (14) 

Quasi-iso. 

60.7  (2.2) 

872  (30) 

The  next  set  of  experiments  consisted  of  constant  amplitude  load-controlled  cyclic  tests 
that  were  conducted  according  to  ASTM  D3479  [4],  Cyclic  testing  was  performed  with  a  constant 
maximum  stress  amplitude,  an  R-ratio  of  0. 1,  and  a  loading  frequency  of  5  Hz.  The  maximum 
applied  stress  amplitude  for  the  90°  specimens  was  varied  between  45%  and  55%  of  the  average 
monotonic  ultimate  stress  of  the  layup  configuration.  The  90°  specimens  failed  by  matrix 
cracking  across  the  width  of  the  specimen.  The  quasi-isotropic  layup  was  tested  with  a  maximum 
applied  stress  amplitude  of  17%  of  the  ultimate  stress  of  the  corresponding  layup. 

Model  Calibration 

The  domain  of  the  microscopic  problem  (i.e.,  RVE)  is  a  unidirectionally  fiber  reinforced 
matrix  as  illustrated  in  Fig.  18.  The  diameter  of  the  fiber  is  set  to  ensure  that  the  volume  fraction  in 
the  RVE  equals  the  experimentally  measured  volume  fraction  of  65.6%.  We  employ  a  4-partition 
reduced  order  model  to  evaluate  the  failure  response  within  the  composite  constituents.  The  matrix 
is  represented  using  three  partitions,  whereas  the  fiber  response  is  idealized  using  a  single 
partition.  The  domains  of  each  partition  within  the  RVE  are  illustrated  in  Fig.  18. 
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Figure  18:  RVE  and  Partition  Structure  for  IM7/977-3 

The  elastic  response  of  the  977-3  resin  is  taken  to  be  isotropic,  with  the  Young’s  modulus 
and  the  Poisson’s  ratio  of  the  material  denoted  by  Em  and  v ,  respectively.  The  IM7  fiber  is 
taken  to  be  transversely  isotropic  with  five  elastic  parameters:  E{ ,  E[ ,  Gl2 ,  v{2 ,  and  vf23 ,  where 
the  1 -direction  is  along  the  fiber  length.  The  Poisson’s  ratios  were  obtained  from  the  literature  ( v\2 
and  v23  from  [12]  and  vm  from  [3 1]).  The  remaining  elastic  parameters  (i.e.,  E[ ,  E[ ,  G{2,  and 

Em  )  were  calibrated  against  the  linear  regions  of  the  stress-strain  curves  recorded  in  the 
monotonic  experiments.  The  calibrated  elastic  parameters  of  the  composite  constituents  and  the 
experimentally  observed  and  simulated  elastic  parameters  of  the  overall  composite  are 
summarized  in  Table  3. 


Table  3.  Elastic  Parameter  Optimization 


Em  [GPa]  E\  [GPa]  E{  [GPa]  G','2  [GPa]  vm  vfl2 

f 

^23 

3.55  263.00  13.00  27.50  0.35  0.32 

0.20 

El  [GPa] 

Ec2  [GPa] 

G,c2  [GPa] 

<2 

Experiment 

158.00 

8.64 

4.66 

0.316 

Model 

158.00 

8.64 

4.66 

0.33 

Superscript  c  indicates  a  composite  material  property. 


The  damage  model  employed  in  this  study  includes  seven  parameters.  Four  of  the  seven 
parameters  (i.e.,  am ,  J3m ,  af,and  /?f)  determine  the  evolution  of  damage  when  subjected  to 
monotonic  loading  conditions,  whereas  the  remaining  three  parameters  (i.e.,  c”,  c”  and  cf) 

determine  the  sensitivity  of  damage  evolution  to  cyclic  loading.  Experiments  conducted  under 
cyclic  tensile  conditions  indicate  that  failure  initiates  within  the  matrix.  Fibers  are  taken  to  be 
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insensitive  to  cyclic  failure  at  the  loading  amplitudes  considered  in  this  work.  The  superscripts  m 
and  f  denote  matrix  and  fiber  phases  respectively.  The  matrix  parameters  are  employed  in 
modeling  failure  in  the  three  matrix  partitions,  whereas  the  fiber  parameters  are  employed  in 
modeling  the  fiber  partition. 


Strain  [mm/mm] 
(a) 


(b) 


Figure  19.  The  Tension  Experiments  are  Compared  to  the  Calibrated  Model  Response 


am  and  a'  are  regularization  parameters  that  control  the  abruptness  of  the  ultimate 
failure  within  the  matrix  and  fiber  phases,  respectively.  The  values  of  am  and  a'  are  chosen  to 
avoid  numerical  difficulties  associated  with  sudden  failure  events,  while  accurately  capturing  the 
characteristics  of  the  stress-strain  response.  /?'"  and  /?'  are  material  parameters  that  control  the 
ultimate  strength  of  the  matrix  and  fiber,  respectively.  The  experimentally  observed  stress-strain 
response  of  the  zero  and  ninety  degree  specimens  are  employed  in  the  calibration  process.  The 
failure  in  the  zero  degree  specimens  is  dominated  by  fiber  failure,  whereas  matrix  cracking 
dominates  failure  in  the  ninety  degree  specimens  subjected  to  monotonic  loading.  J3m  and  /?' 
are  identified  by  minimizing  the  discrepancy  between  the  experimental  and  simulated  stress-strain 
curves.  The  calibrated  model  parameters  are  am  =  0.05 ,  J3m  =  32.0  ,  a r*  =  0.05 ,  and 
/?'  =  340.0 .  Figure  19  illustrates  the  experimentally  observed  and  simulated  stress-strain  curves 
based  on  calibrated  material  parameters  for  zero  degree  and  ninety-degree  specimens.  The  mean 
ultimate  strength  and  strain- at-failure  for  zero  degree  specimens  based  on  experiments  are  2841 
MPa  and  0.0180,  respectively.  The  mean  ultimate  strength  and  strain-at-failure  for  ninety  degree 
specimens  based  on  experiments  are  63  MPa  and  0.00728,  respectively.  The  calibrated  model 
yields  2846  MPa  and  0.0186  for  zero  degree  loadings  and  67  MPa  and  0.00822  for  ninety  degree 
loadings,  which  are  in  close  agreement  with  the  experiments. 

The  cyclic  failure  of  the  matrix  is  characterized  by  the  three  remaining  parameters:  c” , 

c”  and  c” .  The  calibration  of  the  cyclic  loading  sensitivity  parameters  is  conducted  by 
employing  the  stress-life  curves  obtained  from  experiments  in  which  ninety  degree  specimens  are 
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subjected  to  cyclic  loading  in  the  tensile  direction.  In  the  experiments,  the  maximum  amplitude  of 
the  cyclic  loading  was  varied  between  360  MPa  and  520  MPa,  while  keeping  the  R-ratio  constant 
(=0.1).  The  three  parameters  are  calibrated  by  minimizing  the  discrepancy  between  the 
experimental  and  simulated  cycles  to  failure  under  three  different  loading  amplitudes.  A  least 
squares  nonlinear  optimization  algorithm  is  employed  in  the  identification  of  the  optimal 
parameters.  The  calibrated  values  of  c0,  q,  and  c2  are  8.243,  -2.227xl0-2 ,  and  2.192x10  s, 

respectively.  Figure  20  illustrates  the  experimentally  observed  and  simulated  life  curves.  The 
calibrated  model  is  in  close  agreement  with  the  experimentally  observed  mean  stress-life  curve, 
which  is  expressed  in  terms  of  a  power  law  fit.  We  note  that  the  experiments  display  a  substantial 
scatter  around  the  power  law  fit  for  this  material. 


Figure  20:  Experimentally  Observed  and  Simulated  Stress-life  Curves  of  the  Ninety  Degree 

Specimens 


Model  Validation 


The  capabilities  of  the  proposed  multiple  spatio-temporal  modeling  approach  are  validated 
by  comparing  the  predictions  of  the  calibrated  model  with  experiments  conducted  on 
quasi-isotropic  specimens  subjected  to  monotonic  and  cyclic  loading  conditions.  A  quarter  of  the 
specimen  geometry  is  discretized  due  to  symmetry  with  top  eight  plies  explicitly  modeled. 

Figure  2 1  illustrates  the  experimental  and  simulated  stress-strain  response  of  the 
quasi-isotropic  specimens  subjected  to  monotonic  tensile  loading.  The  experimentally-observed 
mean  strength  and  strain-to-failure  are  872  MPa  and  0.0144,  respectively.  The  proposed  multiscale 
model  predictions  of  the  strength  and  strain-to-failure  are  872  MPa  and  0.015 1,  which  are  in 
excellent  agreement  with  each  other.  The  simulations  revealed  progressive  failure  within  the 
matrix  of  individual  off-axis  plies  as  a  function  of  loading.  The  ultimate  failure  is  due  to  fiber 
failure  in  the  0°  plies.  The  predicted  failure  pattern  is  in  close  agreement  with  the  experimental 
observations.  Figure  22  shows  the  results  of  the  acoustic  emission  testing  of  a  monotonically 
loaded  specimen  in  terms  of  hits  as  a  function  of  applied  stress  magnitude.  The  acoustic  emission 
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data  shown  is  indicative  of  the  response  of  a  typical  specimen.  The  distinct  failure  events  predicted 
by  our  simulations  are  indicated  as  well.  The  orientation  of  the  ply  in  which  the  failure  event 
occurs  is  also  shown  in  Figure  22.  All  of  the  failure  events  were  in  the  matrix  with  the  exception  of 
the  final  event.  Since  the  proposed  multiscale  model  is  calibrated  to  the  mean  response  of  all 
specimens,  the  provided  comparison  is  qualitative.  Despite  variations,  the  progressive  nature  of 
the  matrix  damage  accumulation,  as  well  as  the  initiation  of  damage  is  well  captured  by  the 
proposed  multiscale  model. 

We  further  assessed  the  validity  of  the  proposed  model  by  comparing  the  model 
predictions  to  experiments  on  quasi-isotropic  specimens  subjected  to  cyclic  loading  conditions.  A 
sinusoidal  load  with  a  peak  magnitude  of  17%  of  the  mean  ultimate  failure  load  of  the 
quasi-isotropic  specimens  (872  MPa)  and  an  R-ratio  of  0.1  is  applied.  The  objective  of  the 
investigation  is  to  assess  the  capability  of  the  proposed  model  in  capturing  the  distinct  failure 
events  that  occur  within  the  matrix  material  as  a  function  of  the  applied  load  cycles.  In  our 
simulations,  the  first  four  distinct  failure  events  occurred  at  1 1766,  27794,  38817  and  46694 
cycles,  respectively.  Each  of  these  failure  events  was  an  individual  off-axis  ply  losing  load 
carrying  capacity  due  to  a  matrix  crack  extending  the  entire  width  of  the  specimen.  The  first  failure 
event  was  in  a  ninety  degree  ply,  the  second  and  third  events  were  in  fourty-five  degree  plies,  and 
the  fourth  event  was  in  a  minus  fourty-five  degree  ply.  Figure  23a  shows  the  acoustic  emission 
testing  results  in  terms  of  total  hits  as  a  function  of  the  number  of  loading  cycles.  The  failure  events 
predicted  by  the  proposed  multiscale  model  are  indicated  in  Fig.  23a.  Figure  23b  plots  the 
numerical  derivative  of  total  hits  with  respect  to  the  total  number  of  accumulated  load  cycles.  The 
figures  illustrate  that  the  time-to-failure  for  major  recorded  failure  events  coincide  with  those 
predicted  by  the  simulations.  Despite  close  correlation  with  the  acoustic  emission  testing,  the 
predictive  capability  of  the  proposed  model  is  qualitative.  Additional  experimental  investigations 
that  quantitatively  link  the  failure  events  that  occur  during  the  cyclic  loading  are  needed  to  fully 
assess  the  validity  of  the  proposed  approach. 


Figure  21.  Comparison  of  Experimental  and  Predicted  Stress-strain  Curves  of  the 
Quasi-isotropic  Specimens  Subjected  to  Monotonic  Tensile  Loading 
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Total  Hits 


x  104 


Figure  22.  Predicted  Failure  Events  Compared  to  Acoustic  Emission  Data  of  a 
Quasi-isotropic  Specimen  Subjected  to  Monotonic  Tensile  Loading 


Figure  23.  Predicted  Failure  Events  Compared  to  Acoustic  Emission  Data  of  a 
Quasi-isotropic  Specimen  Subjected  to  Cyclic  Loading:  (a)  Hits  vs.  Loading  Cycles;  (b) 
Numerical  Derivative  of  Hits  vs.  Loading  Cycles 
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4.0  RESULTS  AND  DISCUSSION 


This  section  presents  a  two-part  study  on  the  progressive  damage  accumulation  in  carbon 
fiber  reinforced  polymer  (CFRP)  composites.  The  first  part  details  failure  under  monotonic 
loading  conditions.  The  second  part  of  the  study  addresses  fatigue  loadings.  A  multiscale 
computational  homogenization  model  [17]  is  employed  to  numerically  characterize  the 
progressive  damage  mechanisms  of  fiber  fracture,  matrix  cracking,  and  delaminations  as  a 
function  of  loading.  An  experimental  program  using  the  combination  of  AE,  X-ray  radiography 
and  X-ray  CT  techniques  are  employed  to  experimentally  characterize  the  progression  of  damage 
throughout  the  loading  history  and  assess  the  validity  of  the  model  predictions.  A  key  contribution 
is  that  the  sequencing  and  rate  of  failure  at  each  ply  of  laminated  composite  specimens  up  to  the 
sub-microstructure  scale  are  established  based  on  the  combined  experimental  and  computational 
investigation. 

The  remainder  of  this  section  is  organized  as  follows:  Section  4.1  details  the  experimental 
and  computational  investigation  of  CFRP  composites  under  monotonic  loadings  and  Section  4.2 
describes  the  investigation  of  fatigue  loadings. 

4.1  Experimental  and  Computational  Investigation  of  CFRP  Composites  Under 
Monotonic  Loadings 

4.1.1  Macroscale  Numerical  Model 

The  geometry,  finite  element  discretization,  and  boundary  conditions  considered  in  the 
macroscale  specimen  model  are  illustrated  in  Fig.  24.  The  length,  width,  and  thickness  of  the 
numerical  model  were  6  mm,  25  mm,  and  1  mm,  respectively.  The  discretization  of  the  model 
consisted  of  26,560  trilinear  hexahedral  elements.  Each  ply  was  explicitly  modeled  along  the 
thickness  direction  with  16  elements  discretizing  the  thickness  of  the  specimen  layup.  Only  the  top 
half  of  the  specimen  was  discretized  due  to  symmetry  of  the  specimen.  A  small  part  of  the 
specimen  along  the  length  (L=  6  mm)  was  modeled  to  reduce  the  computational  cost  of  the  failure 
simulation.  Periodic  boundary  conditions  were  imposed  along  the  y-direction  to  eliminate 
spurious  boundary  effects  due  to  submodeling.  The  numerical  specimen  was  chosen  long  enough 
to  avoid  the  interaction  of  damage  effects  between  the  top  and  bottom  edges.  The  specimen  was 
loaded  by  increasing  the  average  distance  between  the  top  and  bottom  edges  using  constraints.  The 
magnitude  of  the  applied  stress  was  taken  to  be  the  total  constraint  force  required  to  maintain  the 
specified  average  distance  between  the  specimen  ends  divided  by  the  cross  sectional  area  of  the 
specimen. 
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4.1.2  Microscale  Numerical  Model 

The  symmetric  reduced  order  computational  homogenization  method  presented  in  Section  3.2  is 
used  to  model  the  CFRP  material  response.  In  accordance  with  this  method,  we  describe  our 
choices  for  the  representative  volume  element  and  the  microscale  damage  law  used  to  evaluate  the 
reduced  order  microscale  problem  associated  with  symmetric  reduced  order  computational 
homogenization. 


Figure  25.  The  Unit  Cell  for  IM7/977-3  with  66%  Fiber  Volume  Fraction 


The  unit  cell  of  the  CFRP  composite  material  within  a  single  ply  is  shown  in  Fig.  25.  The 
unit  cell  consists  of  the  unidirectional  fiber  and  the  epoxy  resin.  Consider  the  partitioning  of  the 
unit  cell  domain  into  n  parts  within  which  the  strains  and  damage  are  assumed  to  be  spatially 
constant.  Let  D(u)  be  a  scalar  damage  variable  indicating  the  state  of  damage  within  part  a 
associated  with  the  constitutive  law  in  Eq.  130. 

o(a)  =(l-D(0,))L(a)  :8{a)  (130) 

e (a)  and  <J<a)  are  the  average  strain  and  stress  within  part  a  ,  L(0f)  is  the  tensor  of  elastic 
moduli  of  the  constituent  material  occupying  part  a  ,  and  denotes  the  double  inner  product  of 
two  high  order  tensors.  The  evolution  of  Dia)  as  a  function  of  loading  is  modeled  as: 
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where  t;'"’  is  defined  as: 


(131) 


Z>(B)=®(|£>) 


(“)^  =  —  {u(a,(r)} 


^max  (0  -  max 

0<T<t 


in  which  u(a>  is  the  damage  equivalent  strain  in  part  a  : 


v 


(a) 


lr 


(a)  .  jja)  .  £(a) 


(132) 


(133) 


The  phase  damage  evolution  equation  function  is  modeled  using  a  two-parameter  arctangent  law: 

0(y«) ^  =  arctan(q(Qr)t>(a)  -b(a))  +  arctan(h(0f)) 


71 


+  arctan(h(a)) 


in  which  a(a)  and  b(a)  are  material  parameters  controlling  the  brittleness  of  failure  and  material 
strength,  respectively.  Figure  26  schematically  illustrates  the  effect  of  parameters  a(a)  and  bia) 
on  constituent  material  response. 

The  partitioning  of  the  unit  cell  employed  in  the  present  investigation  is  shown  in  Fig.  27. 
The  partitioning  is  achieved  to  capture  the  three  dominant  failure  modes  of  fiber  fracture, 
transverse  matrix  cracking,  and  delamination.  Parts  1,  2,  and  3  in  Fig.  27  capture  the  liber  fracture, 
transverse  matrix  cracking,  and  delamination,  respectively.  Part  4  is  common  to  transverse  matrix 
cracking  and  delamination.  Introduction  of  this  part  is  an  effective  way  to  treat  intersecting  failure 
paths.  With  this  partitioning  scheme,  the  microscale  reduced  order  model  incorporates  the  relevant 
damage  modes. 
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(a)  (b) 


Figure  26.  Stress-strain  Curves  Produced  by  the  Two  Parameter  Arctangent  Law:  (a)  a(a) 
is  Varied  while  b(a)  is  set  to  Maintain  Constant  Failure  Stress;  (b)  bia)  was  Varied  while 

Maintaining  Constant  a{a) 


Figure  27:  Partitioning  of  the  Unidirectionally  Reinforced  Composite  Unit  Cell 


4.1.3  Calibration  of  the  Model  Parameters 

The  elastic  and  damage  properties  of  the  constituent  materials  (i.e.  fiber  and  matrix)  were 
calibrated  using  experiments  conducted  on  0°  and  90°  unidirectionally  stacked  specimens,  as 
well  as  experimental  data  available  in  the  literature.  In  the  model,  a  uniform  distribution  of  fibers 
was  assumed.  The  variability  seen  in  the  0°  calibration  experiments  is  partially  due  to 
nonuniform  fiber  distribution,  but  as  the  primary  concern  in  this  work  is  tensile  loading,  it  is 
assumed  the  effect  of  nonuniform  fiber  distribution  is  limited. 

The  977-3  resin  was  taken  to  be  isotropic  with  elastic  modulus,  E(m) ,  and  Poisson’s  ratio, 
v(m) .  The  IM7  fiber  was  assumed  to  be  transversely  isotropic  with  elastic  material  properties 
denoted  as  E\i] ,  is®,  G®,  v®  ,  and  v®  .  The  Poisson’s  ratios  of  the  resin  and  fiber  were  set  as 

v(m)  =  0.35  [31],  v,®  =  0.32  ,  and  v®  =  0.20  [12].  The  constituents’  elastic  moduli,  E(m) ,  if®, 
is® ,  and  G® ,  were  calibrated  by  minimizing  the  discrepancy  between  the  composite  elastic 
moduli  of  0°  and  90°  specimens  and  the  simulated  elastic  moduli  of  the  homogenized 
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composite.  The  constituent  moduli  were  determined  as  Eim)  =  3.55  GPa,  T,"1  =  263  GPa, 

Ef  =  13  GPa,  and  G®  =  27.5  GPa  which  were  in  close  agreement  with  previous  investigations 
[31,12], 

The  model  parameters  that  define  damage  accumulation  are  a(m)  and  b{m)  for  the  matrix 
and  a(/),and  b{,)  for  the  fiber.  The  damage  accumulation  parameters  are  calibrated  based  on  the 
set  of  experiments  conducted  on  unidirectional  0°  and  90°  specimens.  The  constituent 
parameters  are  identified  by  minimizing  the  discrepancy  between  experimentally  observed 
stress-strain  response  and  numerical  predictions  of  the  multiscale  model  in  the  least  squares  sense. 
The  calibrated  and  experimentally  observed  stress-strain  response  of  the  0°  and  90°  specimens 
are  displayed  in  Fig.  28.  The  90°  layup  failure  response  was  dominated  by  matrix  failure  whereas 
the  0°  specimens  fail  by  fiber  fracture.  The  matrix  and  fiber  damage  accumulation  properties 
were  therefore  separately  calibrated  based  on  90°  and  0°  response,  respectively.  The  matrix 
parameters  were  calibrated  to  the  maximum  observed  strength  among  the  experimental  scatter 
since  the  failure  originates  at  the  largest  flaw  within  the  resin  along  the  free  edge.  The  fiber 
parameters  were  calibrated  so  that  the  model  response  equals  the  average  experimental  strength 
observed  in  the  0°  specimens.  The  calibrated  model  parameters  were  aUn)  =  0.002 ,  b(m)  =  2.8 , 
a(/)  =  0.05,  and  bU)  =340. 


(a)  (b) 


Figure  28.  The  Stress-strain  Response  of  the  Calibrated  Model  Compared  with 
Experimental  Data:  (a)  Specimens  with  90°  Layup  and  (b)  Specimens  with  0°  Layup 


4.1.4  Results  and  Discussion 

Figure  29  shows  the  cumulative  hit  and  cumulative  energy  as  a  function  of  applied  stress 
amplitude  measured  in  the  AE  testing  at  each  load  increment.  Cumulative  energy  weighs  each 
recorded  hit  based  on  the  magnitude  of  the  strain  energy  released  during  the  damage  event.  When 
the  specimen  is  unloaded  and  reloaded,  the  cumulative  hits  and  energy  remained  relatively  flat 
until  the  past  maximum  loading  magnitude  was  exceeded  indicating  insignificant  cyclic  damage 
accumulation  with  an  exception  between  the  loadings  of  620  MPa  and  710  MPa.  Damage  growth 
initiated  within  the  specimen  indicated  by  an  increase  in  the  AE  hits  around  400  MPa.  Damage 
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within  the  specimen  progressively  accumulated  with  an  increasing  rate  until  the  ultimate  failure  by 
fiber  fracture  in  the  0°  plies.  In  contrast  with  damage  events  at  lower  loading  magnitudes,  the 
acoustic  emissions  at  failure  were  audible  without  any  listening  aides. 


(a) 


Stress  [MPa] 
(b) 


Figure  29.  Loading  Stress  versus  (a)  Cumulative  AE  Hits  and  (b)  Cumulative  AE  Hit  Energy 


While  AE  testing  provides  qualitative  information  about  the  progressive  nature  of  damage 
accumulation,  the  type  and  location  of  failure  associated  with  an  acoustic  hit  is  less  clear.  The 
frequency  and  amplitude  of  recorded  waves  does  provide  some  degree  of  information  on  the 
nature  of  the  failure  event  such  as  fiber  fracture  and  matrix  damage  [14],  but  more  detailed 
information  such  as  damage  in  individual  plies  is  difficult  to  gather  from  AE  measurements  alone. 
X-ray  radiography  and  X-ray  computed  tomography  provides  a  nondestructive  snapshot  of  the 
location  and  type  of  accumulated  damage  within  the  specimen.  An  X-ray  radiograph  of  the  pristine 
specimen  was  taken  before  loading  and  additional  radiographs  were  taken  after  each  loading  (i.e. 
300  MPa,  400  MPa,  620  MPa,  710  MPa,  and  845  MPa).  As  illustrated  in  Figure  30,  the  cracks 
were  visualized  in  light  color  in  the  X-ray  radiographs  due  to  the  presence  of  the  dye  penetrant. 
Since  the  dye-penetrant  could  diffuse  into  the  specimen  through  cracks  originating  at  the  specimen 
edges,  only  edge  cracks  could  be  visualized  in  the  radiographs.  No  substantial  cracks  were  visible 
for  the  first  two  loadings  of  300  MPa  and  400  MPa  other  than  minor  edge  flaws.  At  a  loading  of 
620  MPa,  visible  cracks  appeared  with  orientations  both  perpendicular  and  ±  45°  to  the  length  of 
the  specimen.  Between  620  MPa  and  710  MPa,  the  number  and  length  of  cracks  increased.  Just 
before  ultimate  failure,  a  large  delamination  was  clearly  observed  on  the  lower  left  side  of  the 
specimen  (Fig.  30f).  Smaller  delaminations  were  clearly  observed  on  both  sides  of  the  specimen. 
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Figure  30:  X-ray  Radiographs  after  Loading  to  (a)  0  MPa  (b)  400  MPa  (c)  620  MPa  (d)  845 

MPa 

X-ray  computed  tomography  was  employed  to  obtain  a  3-D  visualization  of  the  extent  and 
mechanisms  of  damage  within  the  specimen.  Figure  31  illustrates  the  3-D  tomographic  image  of 
the  specimen  at  710  MPa  and  845  MPa.  Extensive  45°  and  90°  cracks  are  evident  as  well  as 
delaminations  along  the  length  of  the  specimen  edge.  Figure  32  illustrates  the  layer-by-layer 
damage  profde  observed  using  the  X-ray  computed  tomography  imaging  technique.  The  0°  ply 
shown  in  Fig.  32a  exhibited  some  degree  of  debonding  in  the  fiber  direction.  In  contrast  to  the 
radiography,  the  tomographic  images  are  able  to  capture  damage  zones  away  from  the  edges  that 
are  not  exposed  to  the  dye-penetrant  (Fig.  32a).  The  90°  ply  at  the  center  of  the  specimen 
developed  extensive  matrix  cracking  extending  across  the  specimen’s  width  along  with  rounded 
delaminations  at  the  specimen  edges.  The  45°  ply  (the  top  ply  of  the  specimen)  shown  in  Fig.  32c 
developed  extensive  cracking  across  the  specimen  width  along  with  triangular  delaminations 
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developing  at  the  specimen  edges.  The  large  delamination  shown  in  the  radiograph  of  Fig.  30d 
cannot  be  seen  in  Fig.  32  since  the  tomographic  images  were  taken  over  a  smaller  region  of  the 
specimen  outside  of  the  large-scale  delamination. 


(a)  (b)  (c) 

Figure31.  3D  Tomographic  Images  of  Damage  in  the  Specimen  at  a  Loading  of  (a)  720  MPa 

(b)  845  MPa  (c)  845  MPa 
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Figure  32.  Computed  Tomography  Scans  after  Loading  to  845  MPa:  (a)  0°  Ply  (Ply  11)  (b) 

Central  90°  Ply  (c)  Top  45°  Ply 


The  calibrated  computational  model  described  in  Section  4.1.3  was  employed  to  gain 
further  understanding  of  the  progressive  damage  accumulation  in  the  composite  specimen.  The 
model  provides  a  more  complete  picture  of  the  damage  response  than  the  experiments  alone. 
Figure  33a  shows  the  stress-strain  response  of  the  virtual  specimen  under  mono  tonic  tensile 
loading.  The  predicted  stress-strain  response  of  the  overall  composite  is  displayed  alongside  the 
cumulative  hit  versus  stress  curve  recorded  by  the  AE  system  seen  in  Fig.  33b  (the  hits  were 
summed  over  all  loadings).  The  ultimate  strength  of  the  specimen  predicted  by  the  model  was  855 
MPa  which  was  in  excellent  agreement  with  the  experimentally  observed  strength  of  872  MPa 
with  a  standard  deviation  of  30  MPa.  The  strength  of  the  particular  specimen  probed  by  the  NDI 
techniques  was  846  MPa.  The  ultimate  failure  was  caused  by  fiber  fracture  in  the  0°  plies  in  the 
numerical  investigation. 
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Figure  33.  (a)  Strain  vs.  Stress  for  the  Virtual  Quasi-isotropic  Specimen  (b)  Stress  vs. 
Cumulative  Hits  for  the  Experimental  Specimen 


The  first  major  compliance  change  takes  place  at  approximately  380  MPa  when  matrix 
cracking  initiates  in  the  ±  45°  plies.  The  compliance  change  in  the  virtual  specimen  coincides 
with  the  initiation  of  acoustic  emission  hits  illustrated  in  Fig.  29,  which  occurs  at  approximately 
400  MPa.  The  X-ray  radiograph  (Fig.  30b)  taken  at  400  MPa  displays  insignificant  damage  within 
the  specimen  confirming  the  damage  initiation  prediction  of  the  model.  The  matrix  cracks  that 
initiated  at  the  ±  45°  plies  rapidly  propagate  across  the  length  of  the  specimen.  Figure  34 
demonstrates  the  initiation  and  propagation  of  matrix  cracks  in  the  top  45°  ply.  The  progression 
of  damage  in  the  inner  ±  45°  plies  occurs  less  rapidly  with  cracking  across  the  entire  width  of  the 
specimen  forming  when  the  loading  reached  475  MPa.  The  difference  in  speed  of  damage 
progression  is  due  to  the  confinement  of  the  inner  ±  45°  plies  compared  to  the  top  ply,  which 
retards  crack  growth  when  compared  to  the  top  layer.  Matrix  cracking  within  the  ±  45°  plies  is 
followed  by  the  initiation  of  damage  within  the  90°  plies  between  395  MPa  (when  damage  first 
initiates  in  the  90°  ply)  and  409  MPa  (when  damage  initiates  in  all  90°  plies).  Cracking 
extended  across  the  entire  width  of  the  virtual  specimen  within  the  90°  plies  between  472  MPa 
and  514  MPa.  Figure  35  shows  the  progression  of  cracking  in  the  90°  ply  at  the  middle  of  the 
specimen.  The  90°  cracks  clearly  initiate  from  the  specimen  edges.  Matrix  cracking  in  the  0° 
plies  remains  negligible  until  the  loading  reaches  close  to  the  ultimate  failure  strength  of  the 
specimen.  The  predicted  matrix  cracking  in  the  ±  45°  and  90°  plies  develops  more  rapidly  in 
comparison  to  the  matrix  cracking  in  the  experimental  specimen  as  illustrated  in  Fig.  30b.  This  is 
partly  attributed  to  the  errors  in  the  calibration  of  the  matrix  properties  which  was  conducted  based 
on  90°  unidirectionally  reinforced  specimen  response.  The  failure  in  the  calibration  experiments 
initiated  and  propagated  at  the  most  critical  flaw  (i.e.  the  weakest  link)  within  the  resin.  The 
statistics  generated  by  the  calibration  experiments  therefore  capture  the  lower  end  of  the  strength 
and  ductility  spectrum  of  the  resin  material.  Nevertheless,  the  overall  matrix  cracking  pattern  is  in 
reasonable  agreement  with  the  experimental  observations. 
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Figure  34.  The  Damage  Contours  Corresponding  to  Transverse  Matrix  Cracking  at  the  Top 
Ply  of  the  Specimen  (45° )  at  the  Applied  Stress  Level  of  (a)  388  MPa  and  (b)  394  MPa 


y 


X 


Figure  35.  The  Damage  Contours  Corresponding  to  Transverse  Matrix  Cracking  at  the 
Center  of  the  Specimen  ( 90°  Ply)  at  the  Applied  Stress  Level  of  (a)  412  MPa  and  (b)  470  MPa 

The  initiation  of  delamination  within  the  specimens  occurs  slightly  after  the  initiation  of 
matrix  cracking.  Small  edge  delaminations  initiate  between  385-416  MPa.  The  edge  delaminations 
continue  to  grow  slowly  until  the  loading  reaches  close  to  the  ultimate  strength  of  the  specimen. 
The  rate  of  growth  increases  significantly  as  the  magnitude  of  the  loading  approaches  the  ultimate 
strength.  This  observation  is  in  close  agreement  with  the  high  rate  of  increase  in  the  AE  hits  at  the 
later  stages  of  loading  to  the  progression  of  delaminations  as  shown  in  Fig.  33b.  X-ray  radiographs 
and  tomographs  confirm  this  observation.  Figure  36  illustrates  edge  delamination  propagation  in 
the  7th  ply  around  the  ultimate  strength  of  the  specimen.  The  failure  patterns  predicted  by  the 
model  shown  in  Figs.  34-36  are  in  good  agreement  with  the  patterns  observed  in  the  tomographic 
images  in  Fig.  32. 
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Figure  36.  The  Damage  Contours  Corresponding  to  Delamination  near  the  Middle  of  the 
Specimen  (-45°  Ply)  at  the  Applied  Stress  Level  of  (a)  730  MPa  and  (b)  855  MPa 


4.2  Experimental  and  Computational  Investigation  of  CFRP  Composites  Under  Fatigue 
Loadings 

4.2.1  Macroscale  Numerical  Model 

The  geometry  and  finite  element  discretization  considered  in  the  macroscale  specimen 
model  are  the  same  as  the  model  for  monotonic  loads  discussed  in  Sec.  4.1.1  and  illustrated  in  Fig. 
24.  The  length,  width,  and  thickness  of  the  numerical  model  were  6  mm,  25  mm,  and  1  mm, 
respectively.  The  discretization  of  the  model  consisted  of  26,560  trilinear  hexahedral  elements. 
Each  ply  was  explicitly  modeled  along  the  thickness  direction  with  16  elements  discretizing  the 
thickness  of  the  specimen  layup.  Only  the  top  half  of  the  specimen  was  discretized  due  to 
symmetry  of  the  specimen.  A  small  part  of  the  specimen  along  the  length  (Z=6  mm)  was  modeled 
to  reduce  the  computational  cost  of  the  failure  simulation.  Periodic  boundary  conditions  were 
imposed  along  the  y-direction  to  eliminate  spurious  boundary  effects  due  to  submodeling.  The 
numerical  specimen  was  chosen  long  enough  to  avoid  the  interaction  of  damage  effects  between 
the  top  and  bottom  edges.  The  magnitude  of  the  applied  stress  was  taken  to  be  the  total  constraint 
force  required  to  maintain  the  specified  average  distance  between  the  specimen  ends  divided  by 
the  cross  sectional  area  of  the  specimen.  The  constraint  force  was  oscillated  to  produce  a  fatigue 
loading  with  an  r-ratio  of  0.1  and  a  maximum  stress  amplitude  of  143  MPa. 

4.2.2  Microscale  Numerical  Model 

In  this  study,  we  employ  the  multiple  spatio-temporal  scale  methodology  discussed  in 
section  3  to  evaluate  the  damage  accumulation  and  response  of  the  CFRP  quasi-isotropic 
specimens.  The  unit  cell  of  the  CFRP  composite  material  within  a  single  ply  is  taken  to  be  the  same 
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as  in  Sec.  4.1.2.  Consider  the  partitioning  of  the  unit  cell  domain  into  n  parts  within  which  the 
strains  and  damage  are  assumed  to  be  spatially  constant.  Let  D(a)  be  a  scalar  damage  variable 
indicating  the  state  of  damage  within  part  a  associated  with  the  constitutive  law  in  Eq.  135. 


cr 


(a) 


(l-Z>(a))L(0)  :s(a) 


(135) 


s[a>  and  aw  are  the  average  strain  and  stress  within  part  a,  L<0f)  is  the  tensor  of  elastic 
moduli  of  the  constituent  material  occupying  part  a  ,  and  denotes  the  double  inner  product  of 
two  high  order  tensors.  The  evolution  of  D(a)  as  a  function  of  the  fatigue  loading  is  modeled  as: 


D(a)  =  g‘ 


->(«)  ^  np  d®(p  )  where  0  <  g  —  "  7  <  1 
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where  p  is  the  cyclic  sensitivity  parameter,  (-)+  denotes  the  MacCauley  brackets,  O  is  the 
damage  evolution  law  for  monotonic  loading,  and  uU/  >  is  the  damage  equivalent  strain  defined  as 
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The  monotonic  damage  evolution  law  is  taken  to  be  the  arctangent  law  used  for  the  monotonically 
loaded  specimens. 


(«)\  _  arctan(a(ar)t/a)  -h("))  +  arctan(h(ff}) 
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The  cyclic  sensitivity  parameter  p  is  taken  to  be: 

„(«)  =  Ja)  (a)  (a)  (a)  ,  (a)  >.2 
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where  c0a> ,  c,  ,  and  c2a>  are  parameters  controlling  the  fatigue  response  of  the  material.  u  is 
defined  as: 

^(0  =  max{^(a)(r)}  (140) 

0<T<t 

The  adaptive  macrochronological  time  stepping  method  in  Sec.  3.3.3.  is  utilized  with  a 
A Dmax  =  1% .  The  partitioning  of  the  unit  cell  employed  in  the  present  investigation  is  the  same  as 

that  used  for  the  monotonic  loadings.  With  this  partitioning,  we  can  distinguish  between  transverse 
matrix  cracking,  delamination,  and  fiber  failure. 


4.2.3  Calibration  of  the  Model  Parameters 


The  elastic  parameters  and  the  monotonic  arctangent  damage  law  parameters  da)  and 
h(a)  were  chosen  to  be  equal  to  their  values  from  the  monotonic  loading  calibrations.  The  only 
remaining  parameters  are  c',"’ ,  c\a) ,  and  c\a)  which  control  the  fatigue  response  of  the 

constituent  materials.  A  least  squares  nonlinear  optimization  algorithm  was  used  to  minimize  the 
discrepancy  between  the  mean  stress  life  curve  in  Fig.  37  and  the  calibrated  model  to  determine  the 
fatigue  parameters  for  the  matrix.  The  mean  stress  life  curve  was  detennined  from  tension  tests 
with  fatigue  loadings  on  90°  unidirectional  specimens.  This  is  the  same  mean  stress  life  curve 
from  Sec.  3.3.5  which  should  be  referenced  for  further  information  regarding  its  determination. 
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Figure  37.  The  Calibrated  Fatigue  Response  of  Unidirectional  90°  Specimens  is  Compared 

with  a  Statistical  Curve  Fit  of  the  Experimental  Determined  Fatigue  Response 

The  fatigue  parameters  for  the  matrix  parts  were  determined  to  be  c(0a)  =  8.033  , 
c[a)  =  -6.232 xlO  3,  and  c(2a)  =  -4.35  x  10~6 .  The  fiber  was  assumed  to  accumulate  no  damage 
due  to  fatigue  loadings. 

4.2.4  Results  and  Discussion 

X-ray  radiography  provides  a  nondestructive  image  of  the  type  and  location  of  damage 
accumulating  in  the  specimen.  X-ray  radiographs  of  the  quasi-isotropic  specimens  were  taken  after 
1500,  37500,  72500,  and  100000  loading  cycles.  Figures  38  and  39  show  the  radiographs  from  two 
of  the  three  specimens.  The  two  chosen  specimens  are  representative  of  the  third.  As  illustrated  in 
Figs.  38-39,  cracks  were  visualized  in  light  color  in  the  radiographs  due  to  the  presence  of  dye 
penetrant  applied  to  the  specimen  before  the  radiographs  were  taken.  The  dye-penetrant  diffuses 
into  the  specimen  through  cracks  originating  at  the  specimen  edges,  and  hence,  only  edge  cracks 
are  visualized  in  the  radiographs.  Few  cracks  were  visible  after  the  first  1500  loading  cycles  as 
indicated  in  Figs.  38a  and  39a.  It  is  possible  that  the  cracks  that  are  visible  are  manufacturing  flaws 
as  opposed  to  cracks  that  have  developed  during  the  first  1500  loading  cycles.  However,  it  is 
impossible  to  definitively  draw  this  conclusion  from  the  given  data.  After  37500  cycles,  new 
cracks  have  developed  at  the  edges  of  both  ±45°  and  90°  plies.  By  72500  loading  cycles,  new 
cracks  have  initiated  at  the  specimen  edges  along  with  lengthening  of  some  of  the  cracks  that  were 
visible  at  37500  cycles.  With  few  exceptions  the  cracks  do  not  extend  very  far  across  the  specimen 
with  almost  no  cracks  reaching  half  of  the  specimens’  width.  At  100000  cycles,  the  length  and 
density  of  the  cracks  are  nearly  the  same  as  at  72500  cycles.  No  delaminations  were  seen  in  any  of 
the  specimens  during  the  100000  loading  cycles. 
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Figure  38.  X-ray  Radiographs  of  a  Fatigue  Specimen  after  (a)  1500  Loading  Cycles  (b) 
37500  loading  Cycles  (c)  72500  Loading  Cycles  (d)  100000  Loading  Cycles 
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Figure  39.  X-ray  Radiographs  of  another  Fatigue  Specimen  after  (a)  1500  Loading  Cycles 
(b)  37500  Loading  Cycles  (c)  72500  Loading  Cycles  (d)  100000  Loading  Cycles 

The  calibrated  computational  model  described  in  Section  4.2.3  was  used  to  gain  additional 
understanding  of  damage  accumulation  in  CFRP  quasi-isotropic  specimens  subjected  to  fatigue 
loadings.  We  define  a  transverse  matrix  crack  to  be  present  at  a  material  point  when  the  damage 
scalars  D(1)  and  D(4) ,  associated  with  parts  2  and  4  of  the  partitioned  unit  cell  seen  in  Fig.  27,  are 
greater  than  0.99  .  A  delamination  is  defined  as  when  D(J)  and  D{4)  are  greater  than  0.99  .  A 
fiber  failure  is  indicated  by  D(1)  >  0.99 .  Figs.  40  and  41  show  the  value  of  the  damage  scalar  D(4) 

.  Blue  indicates  a  value  close  to  zero  and  red  indicates  a  value  close  to  1 .  At  every  material  point 
within  this  simulation  D,4)  goes  to  one  after  D(2) ,  and  hence,  in  Figs.  40  and  41  red  indicates  the 
presence  of  transverse  matrix  cracking.  Figure  40  shows  transverse  matrix  cracking  within  the 
third  ply  from  the  top  of  the  specimen  which  has  an  orientation  of  -45°,  and  Fig.  40  shows 
transverse  matrix  cracking  within  the  seventh  ply  from  the  top  of  the  specimen  which  also  has  an 
orientation  of  -  45° .  The  cracking  within  these  plies  was  representative  of  that  in  the  fifth  ply 
which  has  an  orientation  of  45° .  The  top  45°  ply  showed  less  transverse  matrix  cracking  than  the 
internal  plies.  At  1500  loading  cycles,  there  was  no  cracking  within  the  specimen  mirroring  the 
behavior  of  the  experimental  specimens.  After  37500  loading  cycles,  there  was  transverse  matrix 
cracking  within  the  seventh  ply,  but  no  cracking  had  initiated  within  any  other  ply.  By  72500 
cycles,  transverse  matrix  cracking  had  initiated  in  all  ±  45°  plies  with  the  exception  of  the  top  ply. 
At  100000  cycles,  matrix  cracking  had  initiated  in  all  ±  45°  plies.  As  in  the  radiographs,  cracking 
was  limited  to  the  area  in  close  proximity  to  the  edges  of  the  specimen  even  after  100000  loading 
cycles.  At  100000  cycles,  the  simulation  predicted  no  transverse  matrix  cracking  in  the  90°  plies 
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counter  to  the  experimental  observations.  There  was  no  delamination  or  fiber  failure  present  in  the 
simulation  agreeing  with  the  damage  types  seen  in  the  X-ray  radiographs. 


Figure  40.  Transverse  Cracking  in  Third  Ply  (-45°)  of  the  Virtual  Specimen  after  (a)  1500 
Loading  Cycles  (b)  37500  Loading  Cycles  (c)  72500  Loading  Cycles  (d)  100000  Loading 

Cycles 
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Figure  41.  Transverse  Cracking  in  Seventh  Ply  (-45°)  of  the  Virtual  Specimen  after  (a) 
1500  Loading  Cycles  (b)  37500  Loading  Cycles  (c)  72500  Loading  Cycles  (d)  100000  Loading 

Cycles 


The  computational  performance  of  the  adaptive  time  stepping  algorithm  is  shown  in  Fig. 
42.  The  x-axis  is  the  number  of  resolved  loading  cycles  and  the  y-axis  is  the  total  number  of 
loading  cycles.  As  indicated  in  the  figure,  only  966  of  the  100000  loading  cycles  were  resolved. 
This  leads  to  an  average  performance  of  103.5  loading  cycles  per  resolved  cycle.  This  is  a 
substantial  savings  in  computational  effort.  For  the  final  20000  loading  cycles  (cycles 
80000-100000),  the  performance  of  the  algorithm  increased  to  216.2  loading  cycles  per  resolved 
cycle.  Without  the  adaptive  time  stepping  algorithm,  it  would  be  required  to  resolve  all  100000 
cycles  of  a  nonlinear  multiscale  finite  element  model  with  26,560  elements.  This  would  present  a 
difficult  computational  barrier.  The  adaptive  macrochonological  time  stepping  algorithm  allowed 
this  simulation  to  be  run  on  a  single  processor  in  three  weeks. 
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Figure  42.  The  Performance  of  the  Adaptive  Macrochronological  Time  Stepping  Algorithm 

over  100000  Loading  Cycles 
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5.0  CONCLUSION 


This  report  presented  a  multiscale  framework  for  modeling  failure  in  brittle  composites 
subjected  to  monotonic  and  fatigue  loadings.  Numerical  testing  verified  the  framework  which  was 
subsequently  validated  with  experimental  testing  on  carbon  fiber  reinforced  polymers.  More 
detailed  conclusions  regarding  the  multiscale  framework  are  presented  below. 

Section  3.2  presented  a  reduced  order  multiscale  computational  methodology  for  failure 
analysis  of  heterogeneous  materials.  The  proposed  approach  provides  a  novel  model  development 
strategy  for  creating  reduced-order  models  capable  of  efficiently  and  accurately  representing  the 
failure  modes  within  the  microstructure  without  recourse  to  a  detailed  finite  element  model  of  the 
RVE.  A  two-order  modeling  approach  was  devised  to  eliminate  spurious  residual  stresses  upon 
failure  allowing  accurate  stress-redistribution  within  a  macroscopic  component.  The  resulting 
reduced-order  model  possesses  symmetry  allowing  efficient  numerical  evaluation  of  the 
microscale  problem.  The  reduced  order  model  was  verified  against  direct  numerical  simulations. 
The  proposed  model  captures  the  failure  modes  within  the  microstructure  obtaining  good 
accuracy. 

In  Section  3.3,  a  multiscale  computational  framework  for  prediction  of  failure  in  composite 
materials  subjected  to  fatigue  loadings  was  proposed.  The  reduced-order  multiple  spatial  scale 
approach  in  combination  with  the  multiple  temporal  scale  time  stepping  approach  provides  a  high 
level  of  computational  efficiency  without  a  significant  loss  in  accuracy.  This  is  critical  to 
determining  fatigue  life  in  large-scale  composite  structures.  The  experimentally  calibrated  model 
predicted  the  observed  early  life  failure  events  in  carbon-fiber  reinforced  polymer  specimens. 

In  Section  4,  a  comprehensive  experimental/computational  investigation  was  undertaken 
to  determine  the  nature  of  progressive  damage  accumulation  in  CFRP  composites  subjected  to 
monotonic  loading.  Acoustic  emission,  X-ray  radiography,  and  X-ray  computed  tomography 
inspection  methods  obtained  a  clear  picture  of  the  evolution  of  transverse  matrix  cracking  and 
delamination  within  experimentally  tested  composite  specimens  as  a  function  of  the  applied 
loading.  The  multiscale  computational  model  was  employed  to  gain  further  insight  into  the 
interaction  and  sequencing  of  damage  mechanisms  which  were  difficult  to  capture  using  any  of  the 
experimental  techniques.  The  response  mechanisms  captured  by  the  model  predictions  reasonably 
agreed  with  the  experimental  observations. 

An  experimental  and  computational  investigation  was  also  conducted  on  CFRP  composites 
subjected  to  fatigue  loadings  in  Section  4.  X-ray  radiography  visualized  the  evolution  of  transverse 
matrix  cracking  within  the  tested  specimens  as  a  function  of  the  number  of  loading  cycles.  A 
multiple  spatial  and  temporal  scale  model  was  employed  to  simulate  the  accumulation  of  matrix 
damage  within  the  CFRP  specimens.  Considering  the  high  level  of  scatter  in  the  fatigue  calibration 
data,  the  response  of  the  model  showed  reasonable  qualitative  agreement  with  the  experimental 
study.  The  adaptive  time  stepping  algorithm  significantly  reduced  the  effort  required  in  simulating 
fatigue  loadings. 

Within  the  multiple  spatial  scale  framework,  the  challenge  remains  that  the  resulting 
homogenized  macroscale  problem  shows  spurious  mesh  dependency.  The  methodologies 
proposed  in  this  report  eliminate  mesh-dependency  in  the  microscopic  domain,  but  the 
homogenized  macroscale  problem  remains  local.  For  homogeneous  materials,  the  computational 
mechanics  literature  has  extensively  investigated  localization  limiters  that  absolve  spurious  mesh 
dependency  [8],  However,  straightforward  application  of  nonlocal  damage  theory  is  impractical 
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due  to  the  significant  associated  mesh  refinement.  Future  research  will  investigate 
enrichment-based  nonlocal  formulations  that  eliminate  mesh  dependency  without  requiring  a  high 
level  of  mesh  resolution  in  the  macroscale  domain. 

An  aircraft  experiences  a  wide  range  of  temperature  conditions  throughout  its  life  and  even 
within  a  single  flight.  The  properties  of  a  composite’s  constituent  materials  can  vary  significantly 
with  temperature.  These  effects  can  be  taken  into  account  within  the  microscale  domain  where 
constitutive  laws  sensitive  to  temperature  effects  would  replace  the  microscale  constitutive  laws 
presently  utilized  in  the  multiple  spatial  scale  algorithm.  Also,  the  multiple  temporal  scale 
methodology  would  need  extension  to  account  for  temperature  changes  in  macrochronological 
time. 

The  applicability  of  these  methodologies  should  be  extended  from  the  specimen  scale  to 
the  structural  scale.  Several  possibilities  show  promise  in  this  effort  including  the  development  of 
improved  cycle  stepping  techniques  that  require  fewer  resolved  cycles  as  evidenced  by  the 
improved  cycle  stepping  criterion  in  Section  3.3.3.  Also,  the  extension  of  the  proposed  multiscale 
methods  to  shell  elements  would  allow  more  efficient  representation  of  structures  whose  thickness 
is  small  compared  to  its  planar  directions.  In  particular,  this  could  be  utilized  for  modeling 
laminated  composite  aircraft  skin.  Another  path  involves  utilization  of  parallel  computing  since 
the  reduced  order  microscale  problems  at  the  gauss  points  of  the  macroscale  structure  can  be 
solved  simultaneously  within  a  single  time  increment. 

Statistical  variation  is  an  important  consideration  when  trying  to  model  fatigue  failure  in 
composite  materials.  Both  the  constituent  material  parameters  and  the  underlying  microstructural 
geometry  vary  from  specimen  to  specimen  and  even  from  one  part  of  a  specimen  to  another.  The 
scatter  in  the  calibration  experiments  revealed  this  statistical  variability  especially  for  fatigue 
loadings.  In  the  future,  the  multiscale  methodologies  can  be  placed  within  a  probabilistic 
framework  that  rigorously  addresses  and  predicts  the  statistical  variation. 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 


AE:  Acoustic  Emission 

CDM:  Continuum  Damage  Mechanics 

CFRP:  Carbon  Fiber  Reinforced  Polymer 

CHM:  Computational  Homogenization  Method 

CT:  X-ray  Computed  Tomography 

EHM:  Eigendeformation-based  Homogenization  Method 

RVE:  Representative  Volume  Element 

sEHM:  Eigendeformation-based  Reduced  Order  Homogenization  Method  With  Symmetric 
Coefficients 

TFA:  Transformation  Field  Analysis 
UMAT:  User  Material  Subroutine 
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